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Supersonic Diffuser Instability’ 


CHARLES LEE DAILEY# 


University of Southern California 


SUMMARY 


This paper presents the results of an experimental determi 
nation of the nature and cause of an instability that has been found 
to be characteristic of supersonic air inlets operating near maxi 
mum mass flow Steady operation is interrupted by a transient 
phenomenon consisting of a random sequence of individual r¢ 
laxation cveles. Mass flow ecutering the diffuser during steady 
operation is suddenly cut off by a strong interaction between the 
subcritical shock and boundary layer on the surface of the ex 
ternal compression generator which blocks the inlet Air in the 
plenum chamber, stored at high pressure, then ‘blows down" 
until the inlet can restart. The subsequent supercritical flow 
entering the diffuser exceeds the flow rate at the exit, and the 
plenum chamber is recharged to the original condition. A dis 
tinction is drawn between this phenomenon, referred to by the 
author as ‘“‘buzz,"’ and a high frequency wave-type resonance 
noticed at low mass flows and during an individual buzz cvcle 
after the diffuser shock system has been expelled For the 
diffuser described herein, this high-frequency oscillation com 
pares well to the eighth closed-end organ pipe mode of the dif 
fuser at low mass flows and to the ninth mode during the shock 
expelled phase of the buzz cycle. It was also found that buzz 
had a strong adverse effect on combustion stability In most 
cases blowout occurred immediately with the occurrence of the 
first buzz cycle. When burning was maintained during buzz, it 


was found to have no qualitative effect on the buzz cycle 


SYMBOLS 


pressure 

density 

temperature 
ratio of specific heats 
gas constant 

Mach Number 
Reynolds Number 
mass flow rate 

speed of sound 

time 

cross-section area of stream tube 


total internal volume of diffuser model 


Received October 20, 1954 

* The experimental data presented in this paper were obtained 
from several investigations supported by the Power Plants 
Division, U.S. Navy Bureau of Aeronautics 

t Research Associate, Engineering Center 


effective length—ratio of internal volume to plenum 
chamber cross-section area 


a fuel to air ratio for combustion tests 


Subscripts 
() = free-stream conditions 
] entrance to subsonic diffuser 


» 


end of subsonic diffuser 

5) station just downstream of flame arrester 
} end of combustion chamber 

5 sonic station in exit nozzk 

+0) = station at engine exit 

local sonic conditions 

local stagnation conditions 

denotes critical diffuser mass flow 


engine exit 


INTRODUCTION 


T TERE ARE TWO PRINCIPAL requirements ofa super 
sonic diffuser: First it must decelerate the enter 
ing air mass sufficiently to allow steady combustion, 
and second, it must compress the air mass before heat 
addition to achieve good thermal efficiency. 

The former requirement is easily met; it is a simple 
matter to reduce the velocity approaching the burner 
to an arbitrarily low value by appropriate subsonic 
diffusion. 
tion in thrust relative to the burner frontal area. This 


The only penalty is a corresponding reduc 


is offset by an increase in cycle efficiency, but the gain 
is so small that it is of no practical importance. In 
practice, every effort is made to achieve the maximum 
burner velocity compatible with good combustion 
Whatever 
this maximum is, it simply becomes a design require 


efficiency and stable burner performance. 


ment for the diffuser which can easily be met. 

The latter problem is more complicated. In order 
to compress supersonic flow it is necessary to reduce 
the area through which it passes. In order to achieve 
good overall efficiency it is essential to obtain the 
maximum possible supersonic compression. This 1s 
obtained with the least inlet contraction ratio that will 
start. For fixed geometry diffusers, this requires the 
use of a central body to compress the flow before it 
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enters the inlet. Diffusers without central bodies or 
other means of external compression are restricted so 
seriously by the starting problem that they are of no 
practical use at high Mach Numbers. Unfortunately, 
however, the compression produced by the central 
body is accompanied by a cowl drag. Since both drag 
and efficiency increase directly with flow deflection, it 
is necessary to optimize the diffuser design in terms of 
net thrust. Although this is somewhat more trouble 
than designing for a given burner velocity, it can be 
handled in a straightforward manner with existing 
methods. 

An important feature of all supersonic diffusers is 
the existence of the supercritical regime that occurs 
when the inlet volume flow is a maximum for a given 
Mach Number. The “‘critical’’ point is simply the 
point of maximum compression efliciency in the super- 
critical regime and is the condition for maximum net 
thrust of a ram-jet. Operation of a ram-jet at the 
critical point is achieved by increasing heat addition 
in the burner until the shock system within the engine 
is moved upstream to the inlet. The shock losses within 
the diffuser are thus minimized, and the plenum cham- 
ber pressure is a maximum for supercritical operation. 
If further heat addition is attempted after this point 
has been reached, the diffuser enters the “subcritical” 
regime, and either one of two possible situations results. 
Either the flow remains steady, or flow breakdown oc 
curs followed by a nonperiodie fill-up and discharge of 
the plenum chamber known as “‘buzz.”’ 

From a practical standpoint it is obviously important 
to know when this instability might occur, so that 
operation under these conditions can be avoided. It is 
the purpose of the present investigation to determine 
the nature and cause of this instability. 

It should be observed that different kinds of transient 
flow processes may occur in ducts or in internal flow 
systems in general. Depending on the duct geometry 
and flow conditions, these phenomena are sometimes 
even qualitatively different. It would be unrealistic 
to expect a mechanism that accounts for instability in 
one case to apply to all instances of internal flow system 
instability. For example, it would not be expected 
that the well-known phenomena of compressor surge 
and the alternate stalling and unstalling of parallel 
ducts feeding the same plenum chamber could be ex- 
plained by the same mechanism; nor would it be ex- 
pected that axially symmetric supersonic diffusers with 
and without central bodies or scoop diffusers would 
necessarily exhibit the same transient flow character- 
istics. Furthermore, it should not be supposed that 
the mechanism that accounts for instability at a given 
mass flow, for a given diffuser operating at fixed 
Mach Number and Reynolds Number, would be the 
same for all other mass flows. 

It is therefore emphasized that this investigation 
deals principally with the problem of determining the 
cause and nature of the instability that has been ob- 
served to occur in axially symmetric diffusers, with 
central bodies, operating at high suberitical mass 
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flows. 
defines the lower limit of steady subcritical operation 


It is the occurrence of this instability which 


for such a diffuser. 

The first mention of this phenomenon in the literatur. 
appears in an NACA translation of a German report 
by Oswatitsch describing experimental work performed 
by him in 1942. 

He did not discuss the phenomenon in detail but 
merely observed its occurrence. In fact, he dismissed 
the subject with the remark that it was unimportant 
since the subcritical range was impractical anyway 
an observation that is hardly realistic when one realizes 
that a sudden burst of instability might result in the 
failure of a supersonic ram-jet missile or destruction oj 
a high Mach Number turbojet engine. 

This phenomenon was again mentioned by Davidsor 
and Umney in England in 1947, by Ferri and Nucci in 
1948, and is compared to power-plant surge in general 
by Pearce. Stoolman gives a brief discussion of 
Oswatitsch’s experiments and mentions similar results 
that he found in tests of the same size and type of model 
operated over the Mach Number range from 1.8 to 
2.5. However, since Stoolman was concerned prin 
cipally with the possibility of self-excited oscillations 
arising from the intrinsic character of the compression 
system in front of an open duct with low mass flow, 
he did not analyze the high mass flow instability for 
ducts with central bodies. 

In spite of the fact that the occurrence of this phe 
nomenon has been known for more than 10 years, it 
has not yet been adequately explained. The investi 
gation described herein is presented in an effort to give 
a clear physical understanding of this process. 


DISCUSSION OF STEADY OPERATION 


The essential features of steady supersonic diffuser 
performance can be shown quite effectively by the ap- 
proximation of ‘“‘one-dimensional’’ gasdynamics.  Al- 
though the theory is simple and well known, it 1s re- 
produced in detail here for the intuitive feeling, or in- 
sight, which the derivation affords. An understanding 
of the different aspects of steady operation is helpful 
in gaining a physically clear picture of the instability. 

The important features of a supersonic ram-jet engine 
with a central body diffuser are shown in Fig. 1. 


Development of Steady Flow Theory 
The relation between stagnation pressure ratio across 
the diffuser and mass flow is obtained directly from 
the condition of mass conservation. The rate of mass 
flow past any station is 
m = pSM V yRT = (PSM/VT) V y/R 


where y and R are specific heat ratio and gas constant, 
respectively. Mass conservation across the diffuser 


is therefore expressed as follows: 
Py SoMo/V To = P2S2M2/V T>2 


In terms of stagnation conditions, this can be writtet 
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. ryy P rye , , 
So Mo a SoM V Tas Ts 

V1 P/F VT Pelt 
Noting that diffuser flow is essentially adiabatic and 


solving for stagnation pressure ratio, one obtains 


P So Mo V T,,/T. PP 1) 
= ( 
r So Me VT./T, PP 

But since this relation could also be applied between 


upstream and downstream sonic points, it follows that 


P,,/f» = S«,/ ds, (2 


[The condition expressed by Eq. (1) therefore can be 


related to Eq. (2) by a tautology as follows: 
Sx,/ So 


PJ Pe = (3S6/ Se) = 
. Se So 


The area ratio, Sy S», is a dimensionless expression of 
the mass flow entering the diffuser. In practice, it is 
convenient to refer the mass flow to the critical (or 
supercritical) value. Thus one writes the area ratio in 
terms of the critical capture area, Sy,, as follows: 


Sp/Se = (So/So.) (Se-/Ss 


The relative mass flow is now a measure of how close 
the diffuser is to its critical point. The remaining 
factor is a design constant specified for a given engine 
by such requirements as fuel-air ratio, combustion 
velocity, etc. 

The mass conservation condition for a diffuser there 


fore results in the following relation: 


P S /R\S%./S 
P. = - Sa. \ Sa VSe,/ Ss 


Significant Features of Steady Diffuser Performance 


A qualitative picture of diffuser performance as a 
function of relative mass flow can be obtained from 
Eq. (4) and general considerations based on the shock 
configurations in Fig. 2. 

In the supercritical condition (Fig. 2a) and at a con- 
stant free-stream) Mach Number, the mass flow is con- 
stant and the stagnation pressure ratio is fixed by the 
plenum chamber Mach Number. In this condition, 
the stagnation pressure ratio will increase as .\/) de 
creases (resulting from increased heat addition or re 
duced exit area). The maximum value is determined 
by the shock configuration for which the internal shock 
system is as near the inlet as possible. 

If still more heat is added or the exit area is reduced 
further, the first shock of the internal system moves 
outside and intersects the cone shock; the entering 
mass flow is reduced by spillage in the subsonic region 
between the subcritical shock and the inlet. For the 
case shown in Fig. 2b, the shock system efficiency is 
about the same as for the critical mass flow condition. 
The overall diffuser efficiency increases slightly as the 
mass flow decreases because of reduced losses in the 


internal passage. 


I ‘ 
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At low subcritical mass flows, the shock configuration 
shown in Fig. 2c 


this case is that some of the flow entering the diffuser 


occurs. The distinctive feature of 


passes through a single shock of the strong family 
The entropy rise (or stagnation pressure loss) across 
this single shock is greater than that across the two 
shocks inside the intersection. The average stagna 
tion pressure in the plenum chamber therefore de 
creases as the percentage of high entropy air in the inlet 
increases. 

It can be seen on the basis of these qualitative con 
siderations that diffuser performance for steady oper 
ation over a range of mass flows must be shown as 
Fig. 3. 

The dashed line in Fi 
Eq. (4) 


3 is the relation defined by 


for constant exit area (hence constant .\/> 


DISCUSSION OF EXPERIMENTAL DIFFUSER 


CHARACTERISTICS 


No mention has been made up to this point of the 
effect of unsteady operation. All the features shown 
in Fig. 3 have been explained in terms of steady dif 
fuser flow. Although the occurrence in practice of such 
operation is rare, a few cases have been observed. An 
example is shown in Fig. 4a. 

It is interesting to note that the essential features of 
the steady operation performance curve are the same 
as for unsteady operation or buzz. This is shown in 
Fig. 4b where unstable operation commenced when 
the higher entropy air first entered the inlet, so that 
the region of reduced efficiency corresponded to the 
unstable range of mass flows. 

As noted in Fig. 4, the central bodies were slightly 
eccentric relative to the cowls. The fact that entrance 
of the high entropy air coincided with the onset of in 
stability for one case and not for the other might be 
due to the difference in magnitude of this eccentricity. 
It also suggests that this event—-i.e., entrance of high 
entropy flow—may, in some cases, result in instability. 
It is emphasized, however, that the data in Fig. 4a 
show that this condition alone cannot be sufficient for 
instability. The fact that it is not necessary is shown 
by the data in Fig. 5. 

In each case, the stream tube surface passing through 
the intersection of the subcritical shock with the central 
body shock was outside the cowl when instability oc- 
curred. Mach 
Number was too low to allow branching of the strong 


In Figs. 5a and 5c, the cone surface 


subcritical shock near the surface; therefore no entropy 
discontinuities of any kind entered the inlet. These 
cases show that entrance of entropy discontinuities 1s 
not necessary for instability. For the diffuser shown 
in Fig. 5b, the high entropy air following the strong 
shock the cowl when 
instability occurred. 

The linear relation for sonic exit velocity at the 


subcritical was well outside 


throttle area for which high entropy air first entered is 
shown in Fig. 4. Since the slope of the line in Fig. 


tb is much greater than the steady performance curve, 
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2%UPERSONIC DIFFt 


it is apparent that instability could not be caused by 
static divergence of the diffuser following entrance of 
the high entropy air. 

The buzz phenomenon is easily recognized either by 
optical observation or by observation of plenum 
chamber pressure with a system capable of following 
high-speed transients. In wind-tunnel work, the easiest 
way to observe it is by either schlieren or shadowgraph 
observation of the diffuser shock system. The char- 
acteristic features of the oscillation are shown in 
Figs. 6 and 7. 

The optical method not only easily detects the pres- 
ence of buzz but also yields some understanding of the 
phenomenon. It can be seen in Fig. 6, for example, 
that the downstream limit of shock travel during the 
oscillation is the supercritical configuration. Another 
important observation is that the range of shock travel 
does not depend much on the average mass flow. 
In fact, by careful adjustment of the exit throttle, it 
is possible to set the mass flow so close to the transition 
point that several minutes can elapse between indi- 
vidual cycles without affecting the range of shock 
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travel. The interval between oscillations is randon 


the performance of the diffuser during this interval is 
the same as in the steady subcritical range above the 
transition point. By observation of the schlieren 
screen, the period of an isolated cycle can be seen to 
be on the order of 0.2 sec. or less. 
Viewed by itself, this information 
by observation with conventional schlieren or shadow- 


obtained simply 
graph methods—would be little more than suggestive; 
when combined with transient pressure measurements, 
however, a clear picture of the oscillation can be ob- 


tained. 


DESCRIPTION OF MODEL AND EQUIPMENT FOR 
TRANSIENT FLow TEstTs 


The model had a plenum chamber diameter of 5.25 
in. and could be adjusted in length to 96, 120, or 140.4 
in. by inserting cylindrical spools. It was equipped 
with pressure transducers located at three stations on 
the central body and at one station in the plenum 


chamber. 
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Seale drawings of the model, showing the three 
lengths tested, and transducer stations are shown in 
Fig. 8. It had a 20° cone central body, a critical 
capture area at Mach Number 1.95 which was 27.7 
per cent of the plenum chamber area, and a cowl lip 
coordinate angle of 33.8°. 

Cold flow data for this diffuser are shown in Fig. 5b. 
It reached a critical point stagnation pressure ratio of 
0.882 but had a transition point (from steady sub- 
critical to buzz) at 99 per cent of critical mass flow. 

The engine was equipped with a conventional can- 
type burner operating on gasoline fuel. It was in- 
stalled for combustion tests in the 120-in. length engine 
only. The purpose of the combustion tests was to 
determine the influence of buzz on combustion, the 
effect of heat addition on buzz, and transient pressure 
variation during burning and blowout. 

A standard Wiancko system was used for the pressure 
measurements. This consists of a carrier oscillator 
driving an a.c. bridge at 3,000 cycles per sec. The 
two active arms of the bridge consist of coils mounted 
in the transducer. The a.c. current resulting from 
an unbalance of the bridge is demodulated by a full 


wave crystal diode rectifier and used to drive the re. 
cording equipment, which for these tests was a Heiland 
oscillograph with 450 cycles per sec. elements. The 
flat response range of these galvanometers, with damp. 
ing, was only 275 cycles per sec. 

Oscillograph records are presented without reference 
to a pressure scale and are used for a qualitative analysis 
of the transient phenomena. The order of magnitude 
of pressure variation during buzz is shown in Figs, 
12, 13, and 14. 

The wind tunnel had a test section 17 by 20 in. in 
size and was operated nonreturn, without drying or 
temperature control. Air was supplied by a stand-by 
blast furnace compressor at the Kaiser steel mill located 
at Fontana, Calif. This wind tunnel was a part of 
the Fontana Aeronautical Laboratory operated by the 
University of Southern California Engineering Center 
until 1949. This facility has been relocated at 
NAMTC, Pt. Mugu, Calif. It is now known as the 
Aerodynamic Test Division and is still operated, under 
contract from Pt. Mugu, by the Engineering Center. 

The nominal Mach Number of the tunnel was 2.00; 
during the test period it averaged 1.95 with a vari- 
ation of + 0.03 due to humidity effects. The stagna- 
tion temperature of the air stream averaged 200°F. 
with a variation of + 25°. 


DISCUSSION OF EXPERIMENTAL RESULTS 


Typical oscillograms are reproduced in Fig. 9, show- 
ing the effect of average mass flow through the diffuser 
on pressure fluctuations at the 4 transducer stations for 
the intermediate volume. Referring to the figure, 
vertical orientation of the traces down the page corre- 
sponds to downstream orientation of the transducers 
in the model. Vertical timing lines occur at 0.01- 
sec. intervals. 

The important features of the oscillation are clearly 
shown in Fig. 9. At zero mass flow, all stations 
show the presence of a high frequency, periodic os- 
cillation at about 600 cycles per sec. At 18 per cent of 
critical flow, the same high-frequency oscillation 1s 
apparent, but a different phenomenon has now ap- 
peared. A random oscillation with a much higher am- 
plitude and longer period (about 0.035 sec.) now ap- 
pears in addition to the high-frequency motion. At 
41 per cent of critical flow, the long period motion 1s 
dominant; it is more nearly periodic, although still 
random, and the period has increased to about 0.040 
sec. for an isolated cycle. The high-frequency motion 
is still present, although it occurs in bursts at intervals 
that show a consistent relation to the long period cycle. 
Its frequency is now about 700 cycles per sec. At 66 
per cent of critical flow, the long period oscillation has 
acquired a wave form that is characteristic of this 
phenomenon for all higher mass flows between this 
value and the transition point. 

_It appears that the transient operation consists of 
a mixture of two intrinsically different phenomena. 
One is a high-frequency, periodic oscillation that occurs 
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without interruption only at very low flows. The 
other is a relatively long period motion with a fairly 
definite period for each cycle, depending on mass flow, 
but with a random time interval between cycles. Al- 
though the high frequency oscillation is present at all 
mass flows below the transition point, it is not actually 
superimposed on the long period motion. The long 
period motion blocks out the high-frequency oscillation 
over part of its cycle; the interval during which it is 
not blocked occurs consistently at the same phase of 
the long period motion. 

Oscillograms showing the effect of plenum chamber 
volume on the transient pressures just below transition 
are presented in Fig. 10. The second trace from the 
inlet, for the smallest volume case, should be disre- 
garded. The diaphragm of the transducer located at 
that station was broken during the run. 

In addition to the transient pressure measurements, 
schlieren pictures of the inlet shock system were taken 
at a repetition rate on the order of 2,000 to 3,000 
frames per sec. Although the resolution of these pic- 
tures was too poor to permit reproduction in this paper, 
significant observations based on them will be included 
in the following discussion. 

The frequency of the shock motion during the high- 
frequency pressure oscillation at zero mass flow, deter- 
mined by counting frames on the motion picture film, 
was found to be 607 cycles per sec. This agrees with 
the pressure record for zero mass flow in Fig. 9. The 
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short burst of high-frequency oscillation noted above 
for high mass flows can also be seen in the motion pic- 
ture data. The shock pulsation has the same appear- 
ance as in the zero mass flow case, but the frequency has 


increased to 690 cycles per sec. 


Organ Pipe Oscillation at Low Mass Flow 


sefore proceeding with an analysis of the long period 
cycle, it is interesting to see how much can be said 
about the high-frequency motion. 

The physical length of the small volume engine (from 
cowl lip to exit) was 7.83 ft., and the stagnation speed 
of sound was 1,260 ft. per sec. The exit nozzle for the 
high mass flow case consisted of a circular hole in a 1/2- 
in. thick flat plate. The only contraction from the 
plenum chamber diameter of 5.25 in. to the exit diam- 
eter of 2.53 in. was a conical frustum, formed by a 45 
bevel in the 1/2-in. plate with its apex downstream. 
The downstream end of the plenum chamber was there- 
fore essentially a flat plate with a small circular hole in 
the center occupying 23.2 per cent of the base area. 

For the zero mass flow oscillation, the wave length is 
only 2.08 ft. Since the engine length is 7.83 ft., the 
theory of the Helmholtz resonator cannot be used to 
explain this high-frequency motion. It can, however, 
be accounted for in terms of wave propagation. 
the diffuser central body extends only 2 ft. downstream 
from the cowl lip, even the short engine must behave 
At zero mass flow the 
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essentially like an organ pipe. 
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2 304 1.468 5.875 2.510 
2.592 1.520 7.000 2.750 
2 880 1.566 8.000 2 885 
3.168 1.606 8.500 2.925 
3.600 1.650 10.812 3.000 
4.032 1.680 
4.464 1.701 


2° slope to 2.625 








742 JOURNAL OF THE 


| dee ‘ i 
HbAph ery Le ar eue Dod: abaddandi “4 | {4b wer ak: ane bey is 


a 


{ 
vy as vhren* Neha 4ryd 


Liptay hadi doy elt ted phd H 

| } iid yi . AA yh 

YAP, 1 y Tach) “AN J haha: } Pai 
‘hs 4 fh rit hl | t 

TTT TT vane ih ie +] Naat oe) 1) | tee yf beh tea 

(A) SoASoc =O (B) So/Soc = 0.18 


edol,! (My VM At, | PN NA 
aane iF Ly iy % Mide* ity PR Fin AY oh 
bey AL Te Atte TN a ier 
h ral > Or al ~~ ~~, » 
pe s ‘ w / . 

wees us WT NM mn, yer Ww 


Yh SIT ATE Ue tt er 







































































LMT TLE Ie 
i] MNT he i 
| MTT Nal aeereN, 
PAMAT TT ett 
} || f! . 
Soc =0.9! 
0, VN peta henteletnttette 
j | | A ii] | ii] | } i} 1 } 
UU HELL 
T anna eee TM 
| | Midiqi L i 
1] 4 yo Hh Th if an: 
SMPTE LEEAT® LEELA | 















































(6) So/ Soc =0.94 (H) So/Soc = 0.98 


Fic. 9. Effect of mass flow on buzz (effective length, 8.93 ft.). 
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Fic. 10. Effect of plenum chamber volume on buzz cycle. 
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TABLE | 
Comparison of Experimental and Theoretical Frequencies for 
Short Engine 


Experimental Organ Pipe 


Mass Flow Frequency Closed-End Mode Frequency 
0 607 & 604 


During buzz 690 9 685 


exit is closed, and in the high subcritical range is only 
93.2 per cent open. The reflection in the latter case 
is simply an attenuated wave of the same type as for 
the closed-end reflection. The experimental frequencies 
are compared to the eighth and ninth modes of a closed- 
end organ pipe in Table 1. 

The fact that the oscillation is so far from the funda- 
mental mode suggests that the frequency must be 
determined roughly by the driving mechanism. The 
actual frequency is then fixed by the nearest resonant 
mode of the pipe. The forcing function is not under- 
stood. However, a similar high-frequency oscillation 
has been noticed at low mass flows in tests on diffusers 
without central bodies. In both cases, the cowl lead- 
ing edge has been sharp. This suggests vortex shedding 
from the cowl lip as the forcing function, as well as the 
frequency determining factor. It is not clear why the 
frequency should be higher when the exit is partially 
open, although this behavior seems to be characteristic 
of the oscillation. 

Results similar to these have been observed on a small 
diffuser model tested at Mach Number 2. This model 
was only 1.89 ft. long, from the cowl lip to the exit, 
but its transient characteristics were qualitatively 
similar to those of the large engine. 

Because of the small scale, the buzz (long period os- 
cillation) had a much shorter period than that for the 
large engine, but the low mass flow oscillation had es- 
sentially the same frequency. At zero mass flow, the 
high-frequency oscillation occurred at 470 cycles per 
sec. At the test value of the stagnation speed of 
sound (1,160 ft. per sec.) the second mode of a closed- 
end organ pipe 1.89 ft. long is 460 cycles per sec. This 
note was clearly audible. As the mass flow was in- 
creased (below the buzz range), the note suddenly in- 
creased pitch by a little less than an octave, and the 
shape of the shock region in front of the inlet changed 
at the same instant. The actual frequency was not 
measured at the time. However, the third organ pipe 
mode is 767 cycles per sec. (a ratio of 5:3) which agrees 
with the observed change in pitch. Results of these 


two tests are summarized in Table 2. 


TABLE 2 
Comparison of Frequencies for Diffusers of Different Lengths 


Diffuser Experimental Ry, (Based on 
Length, Ft. Mass Flow Frequency Mode Frequency Cowl 
Diameter) 
7.83 0 607 8 604 0.72 X 106 
7.83 >0o 690 9 685 0.72 x 106 
1.89 0 470 2 460 0.74 K 106 
1.89 >0 <940 3 767 0.74 XK 106 





These data indicate there must be some forcing func- 
tion relatively insensitive to frequency in the range from 
400 to 800 cycles per sec., with the actual frequency 
determined by the nearest organ pipe mode. Since 
both frequency and Reynolds Number are essentially 
constant for all cases, it is suggested that vortex shed- 
ding from the sharp cowl lip might be the forcing func- 
tion. It is felt that the experimental evidence justifies 
this as a reasonable hypothesis; verification of the vor- 
tex shedding mechanism through a separate experi- 
mental program would be necessary before it could be 
accepted, of course. 

The main purpose of the present investigation, how- 
ever, is to explain the long period buzz oscillation. 
Although an understanding of the high-frequency os- 
cillation would be both interesting and valuable, it is 
not necessary for the present problem. It is sufficient 
to observe its presence during the phase of the buzz 
cycle in which the strong shock appears in front of the 
inlet. 


Analysis of Buzz 


It can be seen in Fig. 10a that the short burst of high- 
frequency oscillation for the short engine operating at 
high mass flow has a duration of about 0.018 sec. By 
counting frames in the motion picture record it is 
found that the strong shock is in front of the inlet for an 
interval of 0.0185 sec. This shows that the high-fre- 
quency burst, seen in the oscillograms just before the 
inlet pressure trace drops off, corresponds to the inter- 
val during which the strong shock is in front of the inlet. 
Events before and after this interval, shown on both 
the pressure record and schlieren film, can therefore 
be correlated. 

From the film it can be seen that the diffuser shock 
system becomes supercritical immediately after the 
high-frequency oscillation stops. By counting frames, 
it is found that this persists for 0.036 sec., after which 
steady subcritical operation occurs for an arbitrary 
interval. As shown in Fig. 10, the supercritical inter- 
val appears as the low-pressure phase of the cycle for 
the trace nearest the inlet. This is graphically illus- 
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Fic. 15. Effect of fuel-air ratio on burning roughness 


trated by an interval of about 0.012 sec. during which 
the inlet shock is downstream of the first transducer 
station, and the local flow is supersonic. During this 
interval the pressure at this point is a minimum and is 
constant. 

The total fill-up interval—i.e., the period of low pres- 
sure at the first station (Fig. 10a)—is about 0.050 sec. 
Apparently the inlet first becomes subcritical after 
0.036 sec. and reaches equilibrium subcritical oper- 
ation 0.014 sec. later at a slightly higher pressure. The 
buzz cycle has now ended; steady subcritical flow is 
established and persists for an arbitrary time interval 
before another cycle takes place. 

The oscillation is evidently a series of relaxation 
cycles occurring at random intervals between which 
steady, equilibrium, subcritical operation occurs. 

Steady subcritical operation is interrupted in some 
manner, following which the inlet shock moves to the 
tip of the central body, and the high-frequency oscilla- 
tion—characteristic of low mass flow operation—sets 
in. While the inlet shock is oscillating, mass dis- 
charges from the plenum chamber with a corresponding 
drop in pressure. At the end of the high-frequency 
oscillation, the inlet shock moves into the diffuser and 
air enters at the supercritical rate. Because of the 
low plenum chamber pressure, the exit flow rate is less 
than the equilibrium value. The net inflow is there- 
fore positive, and the plenum chamber pressure in- 
creases until the diffuser breaks down again. It has 
been found that, by careful adjustment of the exit 
throttle near the transition point, time intervals of 
several minutes between cycles can be observed. 


NOVEMBER, 1955 


With this qualitative description of the phenomenon 
in mind, a careful analysis of the pressure and schlieren 
data will now be made to determine the cause of the 
instability, the restarting process, and the nature of 
the fill-up phase of a typical buzz cycle. 

It can be seen in the oscillograms that deviation from 
subcritical flow starts first near the inlet. The time 
scale is not large enough to permit accurate determin- 
ation of the transducer that shows the first change in 
pressure. It is clear, however, that the three stations 
on the central body show a pressure drop before the 
plenum chamber station does. The plenum chamber 
Mach Number for this diffuser is 0.117 at the critica] 
point. If this speed were added to the stagnation speed 
of sound (1,260 ft. per sec.), the time for a pressure wave 
to travel from the last central body station to the 
plenum chamber station would be found to be 0.003 
sec., which agrees well with the data. 

Following this initial pressure drop, the high-fre- 
quency oscillation sets in immediately. The time 
interval for the pressure drop cannot be determined 
accurately from these data, but for most cases it ap- 
pears to be on the order of 0.002 sec. This means that 
the inlet flow changes from the steady subcritical con- 
dition to high-frequency oscillation in about 0.002 
sec. This is much less than the time required for par- 
ticle flow through the engine or for wave propagation 
from the inlet to the exit and back. At the plenum 
chamber Mach Number of 0.117, an interval of 0.053 
sec. would be required for a fluid particle to travel from 
the inlet to the exit of the shortest engine. The wave 
propagation time for the same engine would be 0.012 
sec. 

Collapse of the inlet flow—.e., from steady subcritical 
operation to the high-frequency oscillation—has been 
completed in one-third of the time required for an acous- 
tic wave to reach the exit and before a fluid particle 
has had time to enter the plenum chamber. It must 
be concluded from this comparison that conditions 
in the plenum chamber have nothing to do with the 
breakdown of steady flow at the inlet—1.e., the break- 
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Fic. 17. Blowout in steady operation caused by diffuser 
instability. 
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down is so sudden that the way in which it occurs can- 
not be affected by conditions in the plenum chamber. 

The question now arises as to whether conditions in 
the plenum chamber, just before the breakdown of the 
inlet flow, could be the cause of breakdown. Although 
the plenum chamber plays no role in the breakdown 
per se, it is possible that conditions in the plenum 
chamber preceding breakdown may have been sufficient 


to cause the instability. 


Diffuser Separation as a Possible Cause 


Consider the diffuser in a condition of steady subcriti- 
cal operation, and suppose the central body flow were 
to separate suddenly at some point in the subsonic 
diffuser. The central body wake would increase, and 
the core flow velocity in the annular subsonic diffuser 
around the central body would also increase. As a re- 
sult, the pressure at the plenum chamber transducer 
would decrease. If this flow were to mix completely 
in the plenum chamber before reaching the exit nozzle, 
the result would be a decrease in both the local pressure 
and stagnation pressure and a corresponding decrease 
in exit mass flow rate. Mass would therefore begin to 
accumulate in the plenum chamber which would re- 
sult in a compression wave traveling upstream toward 
the inlet. The effect of this wave joining the inlet 
shock would be to move it upstream at the velocity 
required to match the increased pressure downstream 
of the shock. This motion might be unstable, beyond a 
certain amplitude, and result in collapse of the inlet 
flow. 

If this argument were correct, one would expect to 
find a drop in pressure at the plenum chamber station 
as the expansion wave from the separation propagates 
downstream. Suppose, for example, that a separation 
occurred at the last central body transducer station for 
the shortest engine. The pressure would begin to drop 
at the plenum chamber station 0.0030 sec. later. Fluid 
from the separated region would reach the exit 0.046 
sec. later, or 0.043 sec. after the plenum chamber pres- 
sure dropped. An upstream wave would then require 
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Fic. 18. Blowout in steady operation due to rough burning 
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Fic. 19. Blowout during buzz due to buzz. 
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Fic. 20. Blowout during buzz due to rough burning 


roughly 0.007 sec. to reach the inlet. So one would 
expect to find a drop in pressure at the plenum chamber 
station about 0.050 sec. before the inlet pressure 
dropped. For the intermediate and long engines, cor- 
responding values of this time would be 0.067 and 0.078 
sec., respectively. Examination of the oscillograms 
in Figs. 9 and 10 shows that this argument must be 
completely inapplicable—the time intervals are too 
large by a factor of ten to be related in any way to the 
breakdown process. It follows, therefore, that the 
instability is not caused by a transient accumulation 
of mass in the plenum chamber following a loss in dif 
fuser efficiency. 

Since the time interval required for particle flow is 
too long, the question arises whether diffuser flow 
separation could interact on the inlet flow through 
wave propagation in time to cause instability. 

In the diffuser flow separation argument, it was ob- 
served that the exit flow would be reduced when the 
separated flow reached the exit nozzle. Actually, the 
reduction would begin earlier than this. As mentioned 
above, the immediate effect of a central body flow 
separation is a pressure decrease in the core flow, which 
is propagated into the plenum chamber as an expan- 
sion. As the fluid moves downstream from the sepa- 
ration point, a low-pressure region develops between 
that fluid and the separation point on the central 
body—.e., the separated region grows at the particle 
flow rate—which means that expansion waves are 
generated continuously and move into the plenum 
chamber. Consequently, the exit pressure, and hence 
the exit flow rate, begins to drop when the first ex 
pansion wave from the separation reaches the exit 
and drops continuously for some time thereafter. 

It might be argued that mass accumulation in the 
plenum chamber would begin when the first expansion 
wave from the separation reached the exit. In the 
order of increasing engine length, the time intervals 
required for acoustic waves to travel downstream from 
the plenum chamber transducer to the exit and back 
to the inlet are 0.0088, 0.0120, and 0.0153 sec., respec- 
tively. 

If this phenomenon were responsible for breakdown 
of the inlet flow, one would expect to find a smooth 
drooping of the plenum chamber pressure traces be- 
ginning at these time intervals before the sharp drop 
in the central body traces takes place. Furthermore, 
the change in plenum chamber pressure should con- 
tinue until the sharp pressure drop resulting from 
collapse of the inlet flow has had time to be felt at the 





TABLE 3 


Plenum Chamber Pressure Droop Intervals for Three Engine 


Lengths 


Length 
Inlet to Exit, In Ah Ale 44 + Ale 
94 0.0088 0.0038 0.0126 
118 0.0120 0.0038 0.0158 
138.4 0.0153 0.0049 0.0202 


plenum chamber station. In the order of increasing 
engine length these additional time intervals are 
0.0038, 0.0038, and 0.0049 sec., respectively. 

If the separation argument in this form was correct, 
one would expect to find the following situation: The 
plenum chamber pressure trace should begin to drop 
at a time interval, A/,, before the sharp drop represent- 
ing collapse of the inlet flow, and should continue to 
change for an additional time interval, Afo, before the 
sharp drop from the inlet has reached the plenum 
chamber station. In the order of increasing engine 
length, these time intervals are summarized in Table 3. 

Examination of the pressure data in Figs. 9 and 10 
shows these time intervals to be too large by about a 
factor of two. Furthermore, it should be noted that 
the shape of the plenum chamber pressure trace is not 
quite right. Since it is argued that collapse of the 
inlet flow is due to upstream propagation of a pressure 
wave resulting from mass accumulation in the plenum 
chamber, one should expect to find some evidence of 
a pressure rise following the initial drop in the plenum 
chamber transducer trace. Examination of the data 


shows this does not occur. 


Instability Caused by Blocked Inlet Flow 


Another approach to the instability problem can be 
made on the basis of a consideration of inlet flow condi- 
tions alone. In this argument, the plenum chamber 
is not a factor in causing instability, it simply reflects 
the history of the breakdown within the inlet as the 
complete picture, from beginning to end, is propagated 
through the plenum chamber. 

Flow in the inlet of a diffuser is far from the idealized 
concept of a normal shock followed by uniform subsonic 
flow. In reality, it is similar to the shock system in a 
wind-tunnel diffuser, in which transition to subsonic 
flow is accomplished through several shocks. The 
first shock of the train has the characteristic lambda 
form of a strong shock—boundary-layer interaction. 
Near the surface, the flow is separated by an oblique 
shock of the weak family, which forms a Mach reflection 
with a strong shock covering the remainder of the 
entering flow. 

The Mach Number of the flow approaching the 
branch point is roughly 1.59, and since the downstream 
flow is subsonic the intersection will presumably be 
defined by the condition of maximum downstream 
pressure. The maximum pressure ratio across the 
lambda shock occurs when the strong shock is normal 
to the flow. The geometry of this shock for Mach 
Number 1.59 is shown by the solid lines in Fig. 11 for a 
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flow direction of 10° upstream from the intersection 
point. The deflection of 2° across the oblique shocks 
is shown as a separation of the cone surface flow which 
creates a dead water region at the minimum area of the 
inlet arid reduces the effective area available for the 
entering flow. 

A random pressure pulse arriving at this lambda 
shock from downstream would increase the downstream 
pressure above the equilibrium value for the shock, 
Since the strong shock ascending from the Mach re- 
flection is normal, it would have to assume an upstream 
velocity to match the increase in pressure. Calculation 
of the shock intersection problem for the moving shock 
shows a corresponding increase in deflection angle 
across the oblique shocks. An increase of 10 per cent 
in pressure ratio across the normal shock would result 
in an upstream shock velocity that is 15 per cent of the 
local sonic speed and an increase in surface flow sepa- 
ration angle from 2° to 5°. The moving shock is shown 
by the dashed lines in Fig. 11. 

As long as the stream line approaching the cowl lip 
passes through the normal shock, the entering mass 
flow is essentially unaltered. However, as a result of 
both the forward motion and the increased separation 
angle, the effective inlet area available for that mass 
flow decreases as the shock moves forward. 

If a random pressure pulse from downstream has 
sufficient strength to move the lambda _ shock far 
enough forward to block the inlet, mass will accumu- 
late in the inlet, and the subsequent pressure rise will 
As the 
forward velocity of the shock increases, the separation 
At some point, the deflected 


force the shock to diverge to the cone tip. 


angle increases also. 
flow by-passes the entrance entirely, and the entering 
flow is cut off. Spark pictures of the shock during di- 
vergence have shown surface flow deflection angles as 
large as 17°. 

As a result of the abrupt interruption of the flow, 
the moving air column within the diffuser creates a 
sharp pressure drop at the inlet which then propagates 
downstream through the engine. 

The time for divergence of the shock, measured by 
counting frames on the motion picture record, is 0.0022 
sec., which agrees well with the interval of 0.002 sec. 
required for the abrupt drop in inlet pressure in Figs. 
9 and 10. 

This argument would explain why the buzz cycles 
occur more frequently as the mass flow is reduced. 
Since the average position of the inlet shock moves 
forward with decreasing mass flow, a weaker pulse is 
capable of blocking the inlet. Thus a shorter interval 
of steady operation is possible, on the average, before 
such a pulse occurs. 

The question now arises as to why the plenum cham- 
ber pressure droops about 0.005 sec. before the inlet 
pressure begins to drop. An acoustic wave requires 
0.0038 sec. to travel from the inlet to the plenum 
chamber station (for the two shorter engines). This 
means that if the plenum chamber pressure droop were 
caused by something in the inlet, it would have started 
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(0.0088 sec. before the sharp pressure drop at the inlet 
occurred. It appears, therefore, that the inlet blocks 
for about 0.0088 sec. before the entering flow is de- 
flected sufficiently by the advancing shock to by-pass 
the entrance. 

To recapitulate, the instability develops as follows: 
First, the inlet shock is pushed upstream by a random 
pressure pulse from downstream until the inlet blocks. 
Some of the mass that enters the cowl begins to ac- 
cumulate upstream of the point of blockage. The 
pressure in that region rises and results in a forward 
motion of the subcritical shock. This motion is un- 
stable and continues until the flow by-passes the cowl 
lip. After the entering flow has been cut off, the mo- 
mentum of the air column in the diffuser sustains its 
downstream motion for awhile and results in a low- 
pressure region at the inlet. At the same time, the 
mass flow entering the diffuser begins to decrease as 
soon as the inlet blocks and continually decreases as 
the blocking progresses. The pressure downstream 
of the blockage point therefore begins to drop as soon 
as the inlet blocks. It appears from the pressure and 
high-speed schlieren data that the entering flow is cut 
off about 0.009 sec. after the inlet blocks, and about 
0.002 sec. later the shock has advanced to the cone tip 
and the high-frequency oscillation has set in. 

If blockage of the inlet occurred downstream of any 
of the central body transducers, one would expect to 
see a pressure rise at that point coincident with a drop 
in pressure at the nearest downstream station. A good 
illustration of this can be seen in Fig. 9h. The three 
downstream transducers begin to droop at the last 
timing line before the sharp drop occurs, but the first 
transducer shows a rising pressure in this same time 
interval. 

Of the three arguments given above, the last seems 
the most plausible and is in agreement with the high- 
speed schlieren and pressure data. It is believed that 
this explanation is correct. 

Apart from acting as the source of random disturb- 
ances, the subsonic diffuser and plenum chamber have 
nothing to do with the occurrence of instability. The 
breakdown is a local condition in the immediate en- 


trance of the inlet. 


Mechanism of Restarting Supercritical Flow 


The initial pressure drop and subsequent high- 
frequency oscillation have now been identified. Both 
the schlieren and pressure records show an average 
duration of about 0.02 sec. for the high-frequency 
oscillation or about 14 cycles at 700 cycles per sec. 
During this time, the entering flow is cut off, and mass 
discharges from the plenum chamber through both the 
exit and the inlet. 

Referring again to the pressure traces in Figs. 9 and 
10, the next significant event is a large drop in inlet 
pressure caused by restarting of supercritical flow. 
Since the plenum chamber pressure is low when this 
occurs, the inlet flow is supersonic past the first central 
body station and remains supersonic for about 0.012 


sec. As mass accumulates in the plenum chamber, the 
pressure rises, and the shock moves upstream toward 
the inlet. The traces at all stations show steadily rising 
pressure during this fill-up phase of the buzz cycle. 
The pressure increase in the plenum chamber is essen- 
tially due to the increase of stored mass as the diffuser 
approaches equilibrium mass flow. The pressure at 
the first station is the local static pressure at that 
point within the shock region, and it, therefore, increases 
as the shock moves upstream. 

It is interesting to note how well the downstream 
stations reflect the history of the breakdown at the 
inlet—from the initial pressure drop when the inlet 
blocks through the high-frequency oscillation. 

The immediate effect on plenum chamber pressure of 
the inlet restarting appears to be a sharp rise in pres- 
sure. It can be seen in the pressure records (particu- 
larly in Figs. 9f, 9g, and 10b) that this is recorded at 
approximately the same time the abrupt drop in inlet 
pressure occurs which seems to imply that the plenum 
chamber responds to restarting of the inlet instantly 
instead of at the acoustic wave interval of 0.0038 sec. 
The reason for this apparent inconsistency is shown 
partly by the pressure trace at the last central body 
station and partly by the high-speed schlieren data. 

The traces at both the last central body station and 
the plenum chamber station show the same pressure 
jump as the inlet restarts, and the events are separated 
by the proper acoustic interval of 0.003 sec. Since the 
inlet pressure drop, as supersonic flow starts over the 
first central body transducer, coincides with the plenum 
chamber pressure jump, this means that the irlet must 
make an impulsive effort to restart about 0.0038 sec. 
before it becomes completely supercritical. 

Examination of the schlieren motion picture record 
shows this to be correct. The high-frequency oscil- 
lation of the inlet shock continues until supercritical 
flow has started. Transition from the oscillation to 
supercritical operation occurs during the last three 
cycles of the oscillation, which is 0.00435 sec. at 690 
cycles per sec. The first motion of the shock toward 
the inlet occurs at the downstream limit of the excur- 
sion, and entrance of the shock occurs at the same phase 
of the oscillation three cycles later. 

The jump in plenum chamber pressure corresponds 
to the first impulsive entrance of mass into the inlet 
and the inlet pressure drops when supercritical flow 
is established three cycles later. Air in the inlet acts 
like a pump that delivers three impulsive surges to the 
internal air column to get it moving fast enough to 
allow the inlet to become supercritical. 


Calculation of Supercritical Plenum Chamber Fill-Up 


After the inlet has become supercritical, a simple 
approximate calculation of fill-up time can be made, 
based on engine volume and the known rates of flow 
in the inlet and out the exit. In addition to the usual 
assumption of one-dimensional gasdynamics, the ple- 
num chamber pressure is assumed to be uniform and 
equal to the value given by the thermal equation 
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of state for the mass of air in the engine, to be at rest 
filling the entire engine volume, and to be at the up- 
streain stagnation temperature. 

Mass flow, m, through a sonic area, Sx, can be ex- 
pressed in terms of stagnation pressure and tempera- 


ture as follows: 


7+ 1\52- | f, , 
m= | (74 ai Vi [5 ) 
2 R VT, 


The time rate of mass accumulation within the engine 
is the difference between the inlet and exit flow rates. 

a ae 
dm y¥+1\52 l : ’ 
i « si \ x = (SePx — SePr 
dt 2 ev f.. 


= V(dp,,/dt) 
where JV is the total internal volume of the engine. 
By use of the thermal equation of state, the density 
derivative can be written as follows: 


dps, ‘at = (1/RT,.) (dP,,/dt) (4) 


On substitution of Eq. (7) in Eq. (6), one obtains 


= _ (7 : +) 20-0 \ RTs (S,,P 4 SusPa) - 
If an effective length, /, is defined such that 
V = S,/ 
Eq. (8) can be written as 
a(P./P.) _ {r+ a0 Ty On 
l 


(**) Sy ~ S. Py. 
— (9) 
S)oSe \S).S2P 


By intro- 


x 


dt 2 


where subscript, e, refers to the model exit. 


ducing y = 1.400 and rearranging Eq. (9), one obtains 
OUP / Ps) ‘ 1 a, (=) ae 
dt G2 187.47, 


l a (*) Ss | 
O) 
a2)?1\s),S — 


Eq. (10) has a simple exponential solution if the exit 
Mach Number is constant. However, since the ple- 
num chamber pressure is low at the beginning of fill-up, 
the exit Mach Number may be subsonic during the ini- 
tial part of the fill-up period. In the present case, the 
engine extended so far into the wind-tunnel diffuser that 
its exit Mach Number was always subsonic, and Eq. 
(10) was integrated numerically. 

The result of this calculation is compared with a 
typical plenum chamber pressure-transducer curve in 
Fig. 12 for the short engine. If the integration con- 


stant is chosen such that the calculated and experi- 
mental pressures coincide at the beginning of the sharp 
pressure rise at the plenum chamber transducer, it 
can be seen that the calculated pressure rise is signifi- 
cantly less than the experimental curve. The best 
that can be done to bring the calculated and experi- 


mental curves into agreement, within the framework of 
this simple fill-up theory is to assume no exit flow during 
the initial phase of fill-up. This gives the linear pres. 
sure rise shown in Fig. 12 which is still less than the 
data for the first 0.01 sec. 

It appears, from this comparison, that the assump 
tion of uniform density throughout the engine is not 
entirely satisfactory. Immediately after supercritical 
flow has keen re-established, the net rate of mass influx 
isa maximum. The incoming air rushes into the mass 
of air at rest in the engine which, althouch it is com- 
pressed, confines the incoming air to a region near the 
inlet. The initial pressure rise therefore corresponds 
to fill-up in the forward section of the engine rather 
than the entire plenum chamber. The subsequent 
fill-up consists of air entering at a fixed rate but being 
diverted in a nonuniform manner to various regions of 
the plenum chamber as the rather strong compression 
waves shown in the data are propagated back and forth 
in the engine. A more accurate calculation of pres- 
sure variation at the transducer station could be made 
by the method of characteristics, but it is felt that the 
simple fill-up model used here describes the essential 
feature of the phenomenon with sufficient accuracy for 
the present purpose. 

The description of the buzz phenomenon is now 
complete. Steady subcritical operation is interrupted 
by choking of the inlet which causes the subcritical 
shock to move forward to the cone tip and thus de- 
flect the approaching flow outside the cowl lip. As the 
stored mass discharges out both ends of the engine, the 
engine resonates at the ninth mode for a closed-end 
organ pipe of the same length as the engine. As the 
discharge out the inlet decreases, the oscillating shock 
moves toward the inlet and—after three cycles—is 
completely swallowed, and supercritical flow is estab- 
lished. The inlet flow rate is then greater than the 
exit discharge rate, and the plenum chamber fills up 
in an approximately exponential manner until the Shock 
moves out of the inlet and returns to the initial equi- 
librium subcritical operation condition. 

This explanation of the buzz phenomenon is applic- 
able for diffusers with fill-up periods that are long in 
comparison with the period of either the Helmholtz 
oscillation or the fundamental closed-end organ pipe 
oscillation, whichever is greater. If the fill-up period 
is near the period of either of these modes, the relax- 
ation form of the buzz cycle will be distorted by dy- 
namic effects; and if it is less than either of these 
periods, the corresponding mode will appear in place 
of the relaxation cycle even though the downstream 
limit of the shock travel is still supercritical. The 
phenomenon is then similar to buzz with respect to 
the forcing mechanism, but the period is determined by 
the inode (Helmholtz or organ pipe) rather than by 


fill-up time. 


Effect of Combustion 


With the explanation of the buzz phenomenon estab- 
lished, it is interesting to see what effect burning has 
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1 the buzz and also what effect the buzz has on burn 


on t 
ing. 

The engine used gasoline fuel and employed a con- 
Ignition was accomplished 
The pilot 
jet and spark were shut off after ignition of the main 


burner 
operation of the gasoline burner alone and are not in- 


ventional can-type burner. 
with a spark plug and hydrogen pilot jet. 


The data presented, therefore, describe the 


fluenced by re-ignition from the pilot jet. 

A comparison of burner roughness at various fuel—air 
ratios is presented in Fig. 15. Since these tests were 
made with a fixed exit area, the effect of increasing 
fuel flow was to move the internal shocks nearer the 
inlet and decrease the degree of supercritical operation. 
It is possible, therefore, that some of the observed in- 
crease in roughness may have been due to changing 
conditions within the diffuser, as well as to increasing 
burner roughness. If the effect of fuel—air ratio alone 
were desired, it would be necessary to control the exit 
area with changing fuel—air ratio so as to hold constant 
diffuser conditions. 

A comparison of cold flow and hot flow buzz is 
shown in Fig. 16. It is apparent that heat addition 
does not alter the essential features of the buzz cycle. 
On the other hand, it is obvious that buzz affects 
burning, since the airflow through the burner varies 
from nearly zero to the supercritical rate during the 
cycle. 

In most cases, it was found that the extreme vari- 
ations in plenum chamber conditions resulting from 
buzz were too severe to allow the burner to continue 
operating. The only way buzz during burning could be 
obtained was by increasing the fuel—air ratio slowly 
after steady burning was established. In most cases, 
the burner would blow out immediately when buzz 
started. Roughly 20 tries were made before a single 
burning buzz was obtained. 
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Fig. 17 shows a case of blowout caused by the first 
part of a buzz cycle. The droop following blocking of 
the inlet, the sharp pressure drop occurring when the 
inlet is by-passed, the high-frequency oscillation, and 
the entrance of the supercritical shock can all be identi- 
fied. In fact, it 1s interesting to note how well this 
record shows the presence of blocking between the 
first and second transducer stations. The first central 
body trace rises for about 0.01 sec. before the inlet is 
by-passed, and the other two central body traces droop 
during the same interval. The plenum chamber droop 
starts about 0.005 sec. later and has the same duration 
as the central body traces. 

The blowout shown in Fig. 1S was due to rough burn 
ing at supercritical mass flow. The pressure fluctu 
ations at the transducer nearest the inlet are due to 
the shock motion resulting from pressure pulsations in 
the plenum chamber. The blowout appears to occur 
in two explosions; the second is much more violent 
than the first. The upstream propagation of the ex- 
plosion waves is clearly evident in the pressure traces. 
The second explosion was strong enough to blow the 
supercritical shock out of the diffuser, following which 
the high-frequency oscillation occurred for about 
0.004 sec. 

Similar cases of blowout due to kuzz and to rough 
burning are shown in Figs. 19 and 20 where buzz oc- 
In Fig. 19 the burner simply 
When 


supercritical flow was re-established, the flame was out, 


curred during burning. 
quit during the high-frequency oscillation. 


and steady supercritical cold flow ensued. 

In Fig. 20 rough burning appears to occur just after 
fill-up. Apparently the heat release drops off, in- 
creases, decreases, and then increases in the form of an 
explosion. This entire pattern has a definite up 
stream propagation and clearly originates in the ple- 


num chamber. 
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The Stability of Vortex Streets with 
Consideration of the Spread of Vorticity 
of the Individual Vortices 


ULRICH DOMM* 


Hermann-Foettinger-Institute, Technical Unwersity at Berlin-Charlottenburg, Germany 


SUMMARY 

This investigation concerns the stability of vortex streets con 
sisting of two parallel rows of staggered real vortices of identical 
structure whose vorticity spreads with time into the ambient 
fluid. It is found that for the case ot stability the ratio of center- 
line width and pitch of the street is a function of the nondimen- 
sional time r = (4v//?)t, where v = kinematic viscosity, / = pitch 
of the vortex street, and ¢ = time. For vanishing viscosity and 
also for vanishing time, the value calculated by von Karman, 


k =h/l = 0.281, is obtained 


(1) INTRODUCTION 


ly HIS BASIC PAPER on the stability problem of vortex 
streets, von Karman! treated a two-dimensional 
infinitely long vortex street in an ideal fluid in which 
the vortices are arranged on two parallel straight lines. 
He obtained the well-known result that only streets 
with staggered arrangement of vortices are stable and 
that this is the case for one condition only—namely, 
when the quotient of width / to pitch / of the street has 
the value k = h// = 0.281. 

Deviations from this result observed by various in- 
vestigators were explained by Hooker’ by the observa- 
tion that these are due to differences between the 
theoretically assumed and the actual velocity distribu- 
tion, the latter corresponding to that of real vortices. 
If this error is corrected according to Hooker, it appears, 
on the basis of Hooker's experimental data, that e first 
decreases slightly with time and that later it assumes 
values considerably larger than k = 0.281. In a pre- 
vious paper, the author* has calculated the influence of 
the real velocity profiles upon the stability of the single- 
row vortex streets and of the double-row street with 
vortices, arranged oppositely as well as staggered, for 
the special case of perturbation of only one single vor- 
tex. The single-row street was found to be unstable. 
For the double-row vortex street, regions of stability 
were found to exist in the k, 7-plane where k = h / and 
rt = the nondimensional time already referred to. In 
the present paper, the vortex street with staggered 
arrangement will be treated for the case of a perturba- 
tion displacement of all vortices. 

A considerable number of research papers have ap- 
peared concerning vortex streets, and these have been 
reviewed recently by Rosenhead* in a comprehensive 


and stimulating treatise. 
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(2) Bastc ASSUMPTIONS 


According to Hamel’ the spread of a vortex in plane 
real flow may be defined in terms of the kinematic vis- 
cosity v, the time ¢, the circulation at time zero Ty, the 
radius vector y, and the azimuthal velocity component 
€,as follows 


ty 
) (] 


¢ 
arbitrarily small amount from the ideal profile c, = 
For small values of 7 the relation 


For sutliciently large ry and constant ¢, c, deviates an 


To /(2ry). 
cy = [To/(Szvt) ]y 2 


is valid——i.e., there exists a linear increase of c, with 
the radius. The core of the vortex rotates approxi- 
mately like a solid body. 

For the subsequent stability investigation let us use 
the velocity profile, Eq. (1). To establish the equa- 
tion of motion it will be necessary to sum over all indi- 
vidual vortices. Strictly this superposition is not per- 
missible because of the quadratic terms of the Navier- 
Stokes equation. These vanish for the case of the single 
vortex in an infinite field of flow but not in the case of 
disturbed axial symmetry. Let it be assumed that for 
small ratios of core radius to pitch of the street—.e., 
-the center region of the vortex 
In the 


for small values of 7 
may be regarded by approximation as isolated. 
vicinity of the vortex cores, recourse is made to the 
velocity distribution decreasing with | y with conse- 
quent disregard of the quadratic members. The en- 
suing theory may be expected therefore to give results 
at large Reynolds Numbers approximately correct for 
small time intervals and, at small Reynolds Numbers, 


valid also for larger intervals. 


(3) EQUATIONS OF MOTION 


Consider a vortex street, Fig. 1, produced at time 
¢t = Qin the real fluid. This is the model of the vortex 
street treated by von Karman. In addition let it be 
assumed that under the action of viscosity the vorticity 
of all vortices spreads with time into the ambient fluid 
at the same rate for all vortices. The question arises 
whether there exists a function k = k(r) such that small 
impressed perturbations of the vortices remain pre- 
scribably small with respect to the original location of 
the vortices on the parallel straight lines. 
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be (Xm, Vm), those of the lower (x,’, y,’); then there In these equations %) = common translatory velocity 
follows for the radius vector of the mth and mth vortex, of the system. The index m = 0 is not included in the 
relative to the reference vortex (x0, (2) + vo summation; this is incicated by affixing the prime to the 
; summation symbol. 
Noun” = (Xe — Mi — Ly)" + (Pe — Va : 
YOr ' - - Let the perturbations (xm, Ym) and (x,’, y,’) be so 
. feo — [ne + (1/2) — x,"}2 + small in magnitude with respect to the pitch that 
oe ar. linearization of the system [Eqs. (3)] is permissible. 
When introducing the nondimensional time 7 ty /*)t 
The equations of motion for the reference vortex are: this linearization results in: 
ax ¥ —. l = ¢ ‘ I = |” +- (1/2 ] — pb : 
: } By Yo — Ym) + a |; . hs Nad aia ) + 
at STV om m- Sav > iin + (1/2 }? 1 bk 
- 1/2 
2p? e Ps os 2[n + (1/2) ]k 
' Yo — Vn’) 4+ > , : y, &X 
[In + (1/2) ]? +k T Sav “m \y [a + (1/2))]? + Ry 
‘ ‘ n+(1/2)]?4 
; 2[n + (1/2)]k ait ; 
ae ‘ = (Xe — x 
[zn + (1/2)}? + R? T 
| 
dy Po «a, 1 —e€ ""* — 2(t*/r)e To [wn + (1/2)|? —k 
a : (me ~ ta) += 23 b> X 
dr Srv “om m* Srv > > [a + (1/2)]? + Ry 
»b2 ‘ (2 + (1/2) ]?+k?2}/r . 
[a +(1/2)]?+k?2}/r j 2k (é / Mp sa 
ie aye (Yo — Xn) — 2, X 
U[n + (1/2)]? + k? — 25 T Sav “a 
‘ ‘ [ma +-(1/2) ]?+k? r 
2[n + (1/2) ]k = (1/0) 4-88 Je 2[n + (1/2) Jk e ; 
- om | ee — (vo — Vo 
\[n + (1/2)]? + R74? [fn + (1 2)]? + R? T 


1) STABILITY 


Let a, 8, a’, 8’ be functions of the nondimensional time 7. Let the perturbations of the mth vortex of the upper 
row and of the mth vortex of the lower row be of the form 


Xm a ime s a’ i[n+(1/2)]¢ a 
= é R i = ‘ é (.)) 
Vn 6 Vn Pe] 


with 0 < ¢ < 2r. Introducing these formulations into the system |Eqs. (4)] leads to 


a 
a 0 B(r) A(r) C(r) a 
a’ Io B(r) 0 — C(r) —A(r) a’ ( 
= a = (6 
B Srv A(r) C(r) 0 — Bir) |\\ B 
id —C(r) —A(r) —Bi(r) 0 B’ 
with the initial condition 
a(Q) ay 
a’(Q) ao’ . 
= (0a) 
B(Q) Bo 
B'(0) Bo’ 
The equations for a’ and 8’ are obtained by selecting the vortex » = 0 as reference vortex. The functions 


A(r), B(r), C(r); A(r), B(x), C(r) according to Eqs. (4) are defined as follows: 
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1-em | In + (1/2)? — & , 
A(r) = QU’ i-e””)=- —i{i—<¢ , 
2 nl » tin + (1/2)]}? + R}? 
> 2k? e Sian 
e [72 18 2 + } T 
Y / ' l ye 
B(r ae . we ‘ol — 
n | — (1/2 + } \ 
2[n + (1/2)|k ae ‘ 
é 
zm+12)rP +k T 
- [27 + (1/2)]? — k? 
C(r =(, | ee =e a a 
t [2 + (1/2)? + k75* 
2k? é ‘de ; 6b 
t[n-+ 2) le 
+ (1/2)|? + 3 r ) 
e n+(1/2) ]?+k?2}/r e 
A(r) = A(r) + | —>’ (1 em?) | 
T m ¥ j 
a — 2[n + (1/2) ]k : 
B(r) = B(r) = 21 (. : om | Boo a — 
n=0 \)[n + (1/2)]* + R37}? 
2[n + (1/2)Jk eT i tC? Nh: | 
sin | a + - 
[vn + (1/2)]? + R? T 2) * 
e [vn 1/2 T 
C(r) = C(r) —-20 ” allen 
n T 
The mathematical problem is as follows: Find a pi | 1 O 0 i 
system of solutions (a, 8, a’, 8’) of Eqs. (6) and (6a), qi _ ie aan 0 a’ 7 
. of fii ce P F . = } ‘ 
for which (|a|, |8|, |@’|, |8’|) remain prescribably p2 0 - se =I 8 
small as 7 increases steadily. This must apply to all 72 0 0] | p 


a 


y for 0 < ¢ < 2m. To solve, first apply the trans- 


formation 


to the system of differential equations (6), thereby trans- 


muting these to the form 


* 

Pr | B(r) 0 [A(r) — C(r)] 0 pi 

qi = I, 7 0 7 —B(r 0 [A(r) + C(r)] 11 ° 
pe Sry pets) +C(r)] O- B(r) 0 | pe 

qe 0 [A(r) — C(r)] 0 —B(r qe 


It may readily be seen from Eq. (S) that this system 
can be divided into two mutually independent simpler 


systems —namely, 
* 
Pi ly B(r) [A(r) — C(r)]) (pr 
(*") aint Sav UA(r) + C(r)] B(r) f (*") 
(9) 
« 
1 My f — B(r) [A(r) + Clr) ]) fn 
(*) aan Sav \ [A(r) — C(r)] — B(r) f () 
(10) 


with the initial conditions 
pi(0) | — favo + ay’ 
p(0) | \By — av’ 
qi(O) “a a ary’ 
g2(0) Bo + Bo’ 


Since the function B(r) is purely imaginary, Eqs. (9 
and (10) may be simplified further by substituting 


} Ty rR 7. 
S = ) eo Seedy Ban (lla 
fe pe 
le pier 
(a) _ (Ber Eemom ay 
q2 q2 
This calculation yields 
e 
(*) _ I's [ v A(r) — C(r) (?" 
po ’ Sav | A(r) + E(7) 0 pr 
(12 
e 


7 T's [ O A(r) + a (*) 
Je Sav | A(r) — Clr 0 g2 
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Fic. 2. Condition of stability. 


The initial conditions remain invariant when subjected 
to the transformation (lla), (11b). By undoing the 
process of transformation just performed, it is possible 
to trace the path of the individual vortices—namely, 
fora vortex of the upper row of the street 


() 33) 


(28) o* Lave moren + ine] ay 


2 q2 


and for a vortex of the lower row 


( ) - 2 E a) ~ Leer fr Meer 
(0 i N(#) a Lam fi Mmensi(eta)e] (45) 


Since B(r) is purely imaginary, the exponential func- 
tions occurring in Eqs. (14) and (15) cannot exceed 
the value unity. The requirement previously estab- 
lished for the system of solutions (a, a’, 8, 8’ 
namely, that its value remain prescribably small 
may now be applied to the system of solutions (i, 
po, Ji, Go). 

In his previous paper® the author has investigated the 
stability of the solution of a system of differential 
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equations of the type of systems (12) and (13). The 
result of this investigation will be explained for the 


example of Eqs. (12) and (13). The elements of the ma 


} 
trices of Eqs. (12) and (13) are functions of 7, of the 
quotient k = h/l, and of the phase parameter ¢. If 
the matrix elements, say, of the system (12), for a fixed 
value of ¢ are equated to zero, one arrives at certain 
functions k, = k,(r) and k, = k,(r) which derive 
their definition from this process. In the (%, 7)-plane 
both curves, provided they exist, circumscribe a region 
in which the two elements, for any set of (Rk, 7) values 
from this region, are of different sign. In this region 
enclosed by the curves, the solution of the system, in 
the sense discussed above, is stable. The border of this 
region does not belong to the stabie region, except at 
the points that pertain to k, = k,(r) and k, = k,(r) at 
the same time. 

In order to compute the stability range for given 
phase parameters, it is necessary to solve the following 


equations. 


A(r, k, ¢) — Clr, k, ¢) = O 
1(r, k,g) + C(r, k ¢) = 0 

16 
1(7,k, ¢) + C(r,k, ¢) = 0 


Alr, k, ¢ _ C(r, k, ¢ = (J) 


The solution was carried through by numerical and 
graphical methods for the phase parameters ¢: 7 4, 
n/2,39/4, 2, O9/4, 32/2, (4/4. 

It is only necessary to solve the first two equations 
Cl T, k, ¢), 


C(r, k, ¢), by virtue of their definition [Eqs. (6b) ], 


since the functions A(r, k, ¢), A(z, 2, ¢); 


satisfy the relations 


A(r,k, ¢) = <A(r,k, 27 —¢ 
A(r, k, ¢) = i(7,k, 27 —¢ . 
C(r,k, ¢) = —C(r,k, 28 - ¢ ‘ 
C(r,k, ¢) = —C(r,k, 27 -—¢ 


In Fig. 2 the ratio k has been plotted along the ordi- 
nate and the nondimensional time 7 has been plotted 
along the abscissa. Each of the curves of the fam- 
ily of curves plotted has been marked by its phase- 
parameter. First let us consider the stability re- 
gions themselves of the system (12). For the phase 
parameter ¢ = 7 4, the stability range is bounded by 
the curves k,,(7) and k, (7). For r = 0 infinitely 
many solutions are possible. If perturbations with 
phase parameters ¢ — 7m are considered, the possible 
stability range is narrowed steadily. For r = ¢ at 
r = 0, only k = h / = 0.281 is feasible; at r > 0, how- 
ever, there exists a stability range similar to the range 
that results from perturbation of a single vortex. If 
the phase parameter is increased above ¢ = 7, it is 
found that there exist regions of peculiar form, each 
defined by a characteristic point for a certain rt > 0 
which must be common to all stable solutions for that 
particular phase parameter. If a general perturbation 
is considered, composed of perturbations of all per- 
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missible phase parameters, then evidently there re- 
mains only the linear multiplicity of the character- 
istic points as the only possible stable solution. This 
condition has been plotted in Fig. 2 and labeled k = 
k(r). The stable vortex street, in agreement with 
von Karman, for vanishing 7 —-i.e., vanishing viscosity 
or vanishing time—exists at the value k = 0.281. 

The family of curves, Fig. 2, may also be regarded 
as a master diagram of the stability ranges of the sys- 
It is only 
necessary to replace any given index ¢ by 27 — ¢. 
From this it follows in particular that the established 


tem (13), as may be seen from Eqs. (17). 


stability condition satisfies Eqs. (12) and (13) at the 
same time, as required. 


(5) COMMENTS 


The significance of the nondimensional time 7 may 
be understood readily when regarding it as the square 
of the ratio of two lengths. Aside from a constant 
factor, (4vf)'* is equal to the core radius of the real 


vortex, Eq. (1). Although some authors have hinted 


NOVEMBER, 1955 


at the significance of the core radius for the stability 
problem, there exist, outside of Hooker's results, 1, 
experi:nental data. In view of this situation there js 
need for investigations that should furnish information 
whether the kinematic viscosity is adequate for th 
computation of the core radius or else, as may be sus 
pected from experimental investigations,® whether tur 
bulent momentum transport must also be taken into 


account. 
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SYMBOLS 


C conductance ratio = hi//k 
¢, = specific heat at constant pressure 
h = surface coefficient of heat transfer 
k = thermal conductivity 

= half-thickness of wing 
. = Mach Number 
\ position ratio = x// 

complex variable in definition of Laplace transform 
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191] 
SUMMARY INTRODUCTION 
+, pp , ‘ . , : . , P 
he transie : re re distril a vehicle resulting s - ~— , : " — 
rhe transient temperature di tri yution in ; vehicle resulting HE DESIGNER OF AIRCRAFT, mussiles, and other de 
from aerodynamic heating at high supersonic speeds may be a that t hict a i | 
dein rapidly estimated by means of a set of design charts presented — - eS — aren _— x mmniat ™ 
“a herein. These design charts are based on the analytical solution able to obtain reasonable engineering solutions in a 
» to the one-dimensional heat flow in a solid body of constant ma reasonable time to the difficult thermal problems that 
gung terial properties when the boundary conditions are strongly arise in such devices. The thermal problems resulting 
dependent on time Phe surface coefficient of heat transfer is from aerodynamic heating are well known. For ex 
taken to be a linear function of time of flight, and the tempera ha ; : 
ft, p. hapa ; ; ample, a simple calculation shows that for steady flight 
' ture potential is also taken to be a linear function of time of Lape : . 2 
flight. The design charts are presented in such a form that at a speed of 5,500 m.p.h. (or at a Mach Number of 8 
they may be used for cases where the surface coefficient of heat at an altitude of 50,000 ft. the surface temperature of 
trausfer and the temperature potential are much more compli the vehicle will be of the order of 4,000°F. But the 
cated functions of time of flight than used here The method of designer needs to determine rapidly the transient non- 
use of these design charts is illustrated with two examples based i ‘ ‘ . bs ‘ 
ee ee ae uniform temperature distribution from the start of the 
on a previous study of the transient temperature distribution in a ‘ ; stip : ae 
wedge-shaped wing flying at supersonic speeds flight to the steady flight condition or for any specified 
flight program. 
— The general problem of calculation of the transient 


temperature distribution in three dimensions in a de- 
vice moving at supersonic speeds is greatly complicated 
by many factors, such as the geometry, the variation of 
material properties with temperature, the time-depend 
ent boundary conditions, and the state of the art in 
fluid mechanics and heat transfer at high speeds. If 
the problem is restricted to two-dimensional heat flow 
in a solid of simple geometry and with constant materi 


q = rate of heat transfer 

r = integers, 0, 1, 2, properties, numerical methods of integration can be 

{ = temperature employed asin reference 1. If the problem is restricted 

Y = relative time 9 = f ‘ , ‘ , - 
Pastadadamnenden cae further to one-dimensional heat flow in a solid of simple 

v = distance normal to chord of wing ‘ : : 

| : : geometry and with constant material properties, then 

+ = temperature-rise ratio = 1/tp = (t — ty)/(taw — tu) A : z . ir 

this problem with time-dependent boundary conditions 


= defined by Eqs. (36) and (37) 
Z = defined by Eq. (39 


can be solved by analytical means; such a solution is 


a = thermal diffusivity = k/pc, presented here and then is reduced to the form of useful 
eS ; design charts to permit rapid and reasonable estimates 
rT = temperature difference nt : eae : . 

of the transient temperature distribution for many 


= defined by Eq (21 


practical cases of aerodynamic heating. 


i’ = time 

a,b,B,Cy, Ci,Di,D2 = constants In order to demonstrate the utility of these design 
Subccrines charts, the present discussion will use the results ob- 

aw = adiabatic wall condition tained previously for a wedge-shaped wing flying at 

b = initial condition supersonic speeds. These results were obtained by 

p = temperature potential numerical integration of the two-dimensional heat flow 

V = position and also of the simpler one-dimensional heat flow in a 


= wall position 
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solid body with constant material properties. 


STATEMENT OF PROBLEM 


Consider a device moving at variable supersonic 
speeds and assume the local portion under consideration 
is fairly thin, so that one-dimensional heat flow ob- 
tains. As the device moves through the atmosphere it 


rp 
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Fic. 1. Dimensions and airflow pattern of a wedge-shaped wing. 


will rise in temperature, due to aerodynamic heating, 
and heat will be transferred into the device by conduc- 
tion and away from the surface by radiation. Radia- 
tion is neglected here for the purpose of comparison 
with the results of reference 1. The problem is to de- 
termine or estimate rapidly the transient temperature 
distribution in the thin portion of the device for the 
case where the heat-transfer coefficient at a given loca- 
tion varies with time and the adiabatic wall tempera- 
ture varies with time. 

The problem can be illustrated by the wedge-shaped 
wing, shown in Fig. 1, used in reference 1. The prob- 
lem is to determine the transient temperature distribu- 
tion at point e along a line normal to the chord of the 
wing. Fig. 2 shows an enlargement of the wing at 
point e as well as the symbols used. As point e is 
moved to other positions along the wing, a succession 
of temperature distributions normal to the chord is ob- 
tained, and hence an approximate two-dimensional 
transient temperature distribution is produced. 


ANALYSIS 


Consider the general case of transient one-dimensional 
heat flow in a material of constant properties with the 
boundary conditions and the temperature potential 
dependent on time for the simple geometry shown in 
Fig. 2. 

The differential equation is given by 

[O°t(x, 3) /Ox?] — (1/a)[Ot(x, ¥)/Od8] = 0 (1) 


The three boundary conditions are as follows: At 


Atx = 0, 
ot/ox = 0 (3) 


Atx = 1, 
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k(Ot/Ox) = h(tew — te) { 


The surface coefficient of heat transfer and the agj- 
abatic wall temperature are time-dependent quantities, 
given by 


h saat i(d 5 
cas => fold (fj 


For convenience of solution and to save space, di- 
mensionless groups and transformations are used. The 
resulting equations and boundary conditions are then 
solved by means of the Laplace transformation and 
auxiliary methods. 

The following definitions and transformations ar. 
used : 

Temperature rise 


T==i — i r(x, Oo 7 


F4 ce leg — F Sa 


Y = 7/1, (9a 
= (tf — ¢,)/(taw — ty (9b 

Relative time 
X = ad //? (10 

Position ratio 
N=v/ 1] 


Conductance ratio 
C=hli/k (12 
Introducing these definitions, Eqs. (1) to (4) become 


[o27(N,X)/ON?2] — [O7(N,X)/OX] =0 (13 


7r(N,0) = O (14 
07(0,X)/ON = 0 (15 
Or(1,X)/ON = C(X)1,(X (16 


Previous work on the wedge-shaped wing! has shown 
that to a good approximation the variation of surface 
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Fic. 2. One-dimensional heat flow in a wing. 
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coefficient with time may be assumed to be linear and 
that the variation of temperature potential with time 
may also be assumed to be linear. These relations 
will be used merely as first approximations for the 
general one-dimensional heat flow problem. Then 


h=a-+t bd (17 
or C= al/k + (b/* ka)(ad/l?) = Co+ CX (18) 
and tp = BX (19 


where a, b, Co, Ci, and B are constants. 
Hence the boundary conditions given by Eq. (16 


can be rewritten as 
O7(1,X)/ON = (Co + CL.X)BX (20 


The general problem of one-dimensional heat flow, 
defined by Eqs. (1) to (4) with the conditions of time- 
dependence of # and fy. given by Eqs. (5) and (6), has 
been thus reduced to a relatively simple set of equations 
by use of a linear variation of h and 7, with time, given 
by Eqs. (17) and (19). However, it is not necessary 
that A and +, vary linearly with time, since one can 
employ general variations in the form of general poly- 
nomials for the time dependence of # and of r,. The 
resulting general solution, to be published at a later 
date, then uses a double summation for the boundary 
conditions expressed in this analysis by the simple 
quadratic relation of Eq. (20). This general solution 
has also greater applicability than the present analysis 
limited to aerodynamic heating. 

Defining the Laplace transform by 


7(p,N) =| e ?*7(N,X) dX (21) 


and operating on Eq. (13) by this transform, we get the 


following subsidiary equation 


[0°r(p,.N) oN?) — pr(p,N) + 7r(N,0O) = O (22) 


TEMPERATURE DISTRIBUTION 757 


But the initial condition in Eq. (14) permits the drop- 
ping of the last term in Eq. (22), so that it becomes 


[0°z(p,N)/ON?] — pr(p,N) = 0 (23 


The two boundary conditions, given by Eqs. (20) and 


(13) become 


97(p,1) ON = B[(Co,p?) + (2C,/p*)] (24 


07(p,0)/ONV = 0 (25 
Hence the original partial differential equation, the 
initial condition, and the boundary conditions have 
been transformed to a set of three total differential 


equations that can be solved by usual methods. 


9). 
) 


The solution to the subsidiary equation, Eq. (23), is 
given by 
7(p,N) = D, cosh NV p + Dz sinh NY p (26 
07(p,N) ON = Div/p sinh Nop + 
Divi p cosh Nop (27 


Hence 


But from Eq. (25), 

07(p,0) ON = 0 
and we get D. = 0, 
— = Div p sinh Nv p (28 
ON 


From Eqs. (24) and (28) we get 
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D, = [B/(vV/P sinh V/p)] [(Co/p?2) + (2C,/p*)] (29) 


7(p,N) = (B cosh Nv p)/(sinh Vp) X 
[(Co/p’*) + (2C,/p"’*)] (30) 


Hence 


In order to transform Eq. (30) back to the N-X 


plane, it was found necessary to change the hyperbolic 


functions in Eq. (30) by the binominal expansion into an 


infinite series as follows: 


> 


lev o(1 — e-2 


— 


The desired inverse Laplace transform given on 
page 380 of Conduction of Heat in Solids by Carslaw 


cosh Nv p/sinh Vp = (eVV?+ e-NvVP) + 
\ p)] 
(31) 


= > le-Vve e+! V) + e-V > (2 


and Jaeger is now written in terms of p, Y, and the 


parameter x, as follows: 


jr l 
l2y xs 


ila, 


tp! + n/2§ 


£ = (4X)”* 2” erfe ; (32) 


For the rth term of the infinite series which is formed 
by combining Eq. (31), with the first term in the first 
set of parenthesis in Eq. (30), we get 


e-v?P (2r + 1 V) e-v?e (2r +1+N | 
£' {BC E+ 372 + 1 + 3/2 It ~ 
P P 


(2r + 1 — N) 
2/X 
(27+ 1+ wt 
24/ x 


BO ; 4x)" E erfe 


t erfe (33) 
For the rth term of the infinite series which is formed by 
combining Eq. (31) with the second term in the first 
set of parenthesis in Eq. (30), we get 


f e—~V Pe (2r +1 V) e-veCr+t+N \ 
£-' <2BC, il + 5/2 ~ + 1 + 5/2 7 


" - p 
- eet. . tart i — NN) 
2BC, 44 y”"" | ® erfe 2./X 
2V/. 
 (Q@r+i4 ay 
a” eric ae ; 34) 
2X f 


The complete solution, transformed to the -X 
plane, is given by the following infinite series: 


t(N,X) = BY [(Co(4X)*7 (a8 erfe s + 7% erfe 2’) + 
r=0 
2C,(4X)'°(@ erfe s + ® erfe s’)} (35) 


where 2=(2r + 1 — N)/2VX 
2’ = (2F + 1+ N)/2VX (37) 


The analytical solution represented by Eq. (35) can 
be verified by showing that it satisfies simultaneously 
the original differential Eq. (13) and the three bound- 
ary conditions [Eqs. (14), (15), and (16)]. The solu- 
tion can also be used to obtain the spatial and time 
temperature gradients; these will not be given here but 
in a later publication discussing the generalized solu- 
tion for any variation of boundary conditions with time. 
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The complete solution given by Eq. (35) can be 
placed in a more useful form by eliminating the con 
stant B of Eq. (19) since it is not zero and introducing 
the temperature-rise ratio of Eq. (9), as follows: 


y (ft = to) {a — by) 


II 


C22 X & (8 erfe s + 7 erfe 2’) + 


’ 0 
9/90\87 v3/2 “s se, 5 ag 9 
2(2)5 CX > (® erfes + Perfes’) (38 
r=0 


Let us next define, for brevity, a function Z,,, which 
includes the infinite sum of the repeated integrals of the 
error functions given in Eq. (38), 


, m — | > pm oe ee 
Ln = 2m ( : ) v/ X =: (@”" erfes +7” erfes’) (39 
o r=0 


where m is an integer indicating the order of the re 
peated integral of the error function. 

The final result for the temperature distribution can 
be written as 


Y = CoZ3 + CXZ; +) 


where the constants Cp and C are defined in Eq. (15). 

The transient temperature distribution, given by 
Eq. (40), can be computed in a straightforward way if 
the various Z,, functions are computed and arranged 
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DESIGN CHARTS FOR TRANSIENT 


in chart form. The next section describes these func- 
tions and illustrates their use for problems of aerody- 


namic heating. 


DESIGN CHARTS FOR TRANSIENT TEMPERATURE 
DISTRIBUTION WITH TIME-DEPENDENT BOUNDARY 
CONDITIONS 


Tbe major objective of this paper is to present a set 
of useful design charts to permit the designer to predict 
the transient temperature distribution for one-dimen- 
sional heat flow in a thin body moving at supersonic 
speeds under conditions of aerodynamic heating. 
These design charts are based on the analytical solution 
given by Eq. (40), where the temperature-rise ratio, Y, 
isa function of the relative time, X, the position ratio, 
\, and certain arbitrary constants, dependent on the 
flight conditions. Since Z,, is a function only of rela- 
tive time, XY, and position ratio, .V, it is simple to con- 
struct charts, each for a different value of NV showing 
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many values of Z for a large range of values of X. Such 
charts have different coordinates than those used in the 
well known Gurney-Lurie charts and the Hottel charts, 
constructed for the special case of boundary conditions 
that are independent of time. 

Six design charts have been constructed for values of 
the position ratio, NV, equal to 1, 0.8, 0.6, 0.4, 0.2, and 0, 
shown in Figs. 3, 4, 5, 6, 7, and 8, respectively. Each 
of these six design charts contains the Z,, functions, 
Z;, Z;, Z;, Z:, Zs, and Z,,, so that the charts may be 
used to solve problems that are more general than 
those considered in this paper. Various attempts 
were made to reduce the several lines for different 
values of Z,, to a common line, but in all cases tried a 
considerable loss in accuracy was introduced. The 
charts were computed with the aid of a recent tabula- 
tion of repeated integrals of the error function.” 

The six design charts are presented in log-log form 
for several reasons. First the state of the art of heat- 
transfer coefficients for supersonic flow does not permit 
selection of this important parameter with an accuracy 
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of better than 20 per cent. Second, the properties of 
the material are assumed to be independent of tempera- 
ture, but actually this factor introduces additional error 
at the higher temperatures. Third, the heat flow is 
not truly one-dimensional in a real case. Hence, it was 
decided to select scales for Y and X on these charts 
which would permit engineering estimates of fair ac- 


curacy. 


Methods of Use of Design Charts 

The design charts are used most easily by evaluating 
first the surface temperature of the device at the de- 
sired instant of time and then the internal temperatures 
by the following scheme. 

At the surface, V = 1, and so 


VY, = 11/tp = 11/(Taw — 71) (41) 
If Eq. (41) is solved for 71, we get 
7 = TawV1/(1 + Vi) (42) 
For any internal position ratio, V, we have 
Vy = ty / tp = tw/ (tae — Ti) (43) 
Considering Eqs. (42) and (43), we get 
™N = Taw Yn /(1 + V1) (44) 


From the conditions of the types of problems con- 
sidered here, X is computed first, the six charts entered 
at this value of XY to produce the values of the desired 
functions of Z,,. Then the value of Y;, at the surface, 
is computed from the final solution given by Eq. (40), 
or by variations of this solution. The surface temper- 
ature rise is computed from Eq. (42), once the adiabatic 
temperature rise, T,., has been found from the condi- 
tions of the problem. The individual values of Y,, are 
then computed from the known values of Z,, for each 
value of the position ratio, N. Finally, Eq. (44) is 
used to compute the temperature rise ry at each value 
of N. Knowing the initial temperature, /,, the actual 
temperature distribution is then found. 

The method of using the design charts is illustrated 
by two detailed numerical examples. In order to save 
space these examples are selected from the previous 
study! for the transient temperature distribution in a 
wedge-shaped wing flying at supersonic speeds. In 
this study a wedge-shaped steel wing, with assumed 
constant properties, is initially flying steadily at a 
Mach Number of 1.4 at a constant altitude. It is 
suddenly subjected to a constant acceleration until it 
reaches a Mach Number of 6 and then is flown steadily 
at a Mach Number of 6. The two-dimensional heat 
flow in this wing was solved first in this study to obtain 
the transient temperature distribution in the wing. 
It was also shown that the two-dimensional heat-flow 
problem could be reduced to the solution of a series of 
one-dimensional heat flows at various positions along 
the chord of this wing. Both types of solutions were 
based on integration of the heat-flow equations by 
laborious numerical methods. The present analytical 
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solution has been sought as a means of reducing sig- 
nificantly the labor of computations to get reasonable 
answers in a reasonable time. 

The reader is referred to the previous study for the 
values of the properties of steel, for the variation of sur- 
face coefficient of heat transfer with position and with 
time or Mach Number, and for other details such as 
recovery factors used and variation of adiabatic wall 
temperatures with speed. For this reason only a mini- 
mum amount of repetition will be used in the examples 
taken from this study to illustrate the present analytical 


solutions. 
EXAMPLES ILLUSTRATING USE OF DESIGN CHARTS 


Example 1 

Consider the wedge-shaped wing of the previous 
study! shown in Fig. 1 flying steadily at a Mach Num- 
ber of 1.4 at an altitude of 50,000 ft.; suddenly it is 
accelerated to a Mach Number of 6 at a constant ac- 
celeration of 2 Mach Numbers per min. and then flown 
steadily at a Mach Number of 6. The half-thickness 
of the wing at 2 ft. from the leading edge is 0.10 ft. 
The problem is to determine the transient temperature 
distribution normal to the chord at 2 ft. from the lead- 
ing edge as the Mach Number changes from 1.4 to 6. 
For purposes of illustration, only the temperature dis- 
tribution when the wing first reaches a Mach Number 
of 6 will be computed here, since the analytical solu- 
tion can be applied immediately to solve the problem 
for any instant of time without going through the en- 
tire history of solutions up to that time, as required in 
the laborious numerical integrations. 


Solution by the Design Charts 

From the previous study, the variation of the surface 
coefficient of heat transfer at 2 ft. from the leading edge 
is given by the approximate expression, 

hb - S57 15 
But VW AT, + dd b 
where d is the acceleration. Hence, combining Eqs 
(a) and (b) 

h = 20 + 586 c 
where # is in seconds, and Eq. (c) is of the form of Eq 
(17). When the properties of steel are inserted in Eq 
(18), together with Eq. (c), we get 
C = al/k + (b73/ka)X 


= 20 X 0.1/22 + 


[(5/6)(0.1 X 0.1 X 0.1 & 3600) /(22 XK 0.384)] X 
= 0.09091 + 0.3551 X d) 
where Cy = 0.09091 and C; = 0.3551. The time of 


flight is (6.0 — 1.4)/2 = 2.3 min. Hence 
X = ad/|? = (0.384 X 2.3)/(0.01 & 60 


= 1.472 
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0.067 Approx. M.P.( °F) 1250 2350 2750 3300 
0 k (Btu/ft hr °P) 70 92 22 8. 
10 0.4072 0 2 0.2746 2412 0.22 fF (1vm/rt?) 172 114 486 276 
d) + (ee - 462.0 290.9 24 zz7 _(Btu/lbm °F) 0.23 0.52 0.118 0.13 
° 7 ] z 462 4] a (ft</nr) 1.76 1.5 0.384 0.246 
e of l 1 
' 4ye sae 2 a 0.1183 0.1022 1.726 6.762 
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TABLE 3 
Temperature Distribution At Instant Mach © is Reached in 


Wing of Example 2 








Aluminum Beryllium Steel Titanium 
x 2.999 2.641 0.654 0.4174 
N21 1.320 1.200 0.494 0.390 
0.8 1.148 1.026 .330 .231 
0.6 0.995 0.874 .221 Be ey 
Ze 0.4 .915 .800 152 . 0808 
0.2 .853 .734 113 051 
(@) Buy 728 .0975 040 
N=1 0.4684 0.3240 0.5578 2 
0.8 -4074 .2768 3726 0.65 
Y= 0.6 .3531 .2360 2498 
C,XZe 0.4 . 3247 .2160 .1716 .2 
1 0.2 .3027 .1982 ‘1276 "1454 
0 .2995 .1966 .1101 .1143 
N=1 685.5 525.9 769.5 1126.1 
0.8 596.2 44g, 514.0 667 .0 
0.6 516.8 383.0 344 2 395 .€ 
* 0.4 475.2 350.6 236 .8 33.3 
(FP) 0.2 443.0 321.7 176.0 148.7 
0 438 .3 319.1 151.9 116.9 
N=1 757 598 841 1198 
0.8 668 521 586 739 
t(°P) 0.6 588 455 416 467 
0.4 547 422 308 305 
0.2 515 393 248 220 
0 510 391 224 189 





The desired solution is of the form of Eq. (40). Hence 
Y = (4; + CXZ; (40) 
so that Z; and Z; are found from each of the six design 
charts, Figs. 3-S, for the above value of XY, and then the 
value of Y is obtained from Eq. (40) at each of the six 
position ratios in the wing. The values of Y are then 
combined as in Eqs. (42) and (44), to obtain the tem- 
perature rise, 
7; = Taw V1 (1 + Y}) (42) 
TY = Taw YN (1 + Y;) (44) 
from the adiabatic temperature potential, for a recov- 
ery factor of 0.9, as used in the previous study! 
Taw Meek = 1d 


2,149°F. 


II 


These calculations are summarized in Table 1. 

The results of the numerical integration for the same 
problem, on a one-dimensional heat-flow basis, are 
shown in the bottom row of Table 1, labeled ¢’. The 
temperatures calculated with the aid of the design 
charts agree with the previously calculated tempera- 
tures, ¢’, within the accuracy of the reading of the pres- 
ent design charts and within the accuracy of the finite- 
difference process used in the numerical integration 
with the relatively coarse network. For engineering 
purposes, both sets of temperatures represent excellent 
answers to the problem, but those based on the charts 
required only a small fraction of the time used to get 
the results based on the numerical integrations of the 
one-dimensional heat-flow equations. 

The values of the Z functions shown in Table 1 were 
read from the original master charts with a maximum 
error of 2 or 3 digits in the third significant figure when 
the value of Z is less than unity. However, the values 
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of the Z functions, read from Figs. 3 to 8, will have 
maximum reading errors of the order of five to eight 
digits in the third significant figure of Z. Where 
greater reading accuracy is warranted, tables of the Z 
functions with proper interpolation rules can be pre 
pared. 

It is evident that the analytical solution, based on the 
one-dimensional heat flow, can be used to predict with 
excellent accuracy the transient temperature distribu- 
tion in the wing which results from two-dimensional 
heat flow, except near the leading and trailing edges 
This result is based on the work done in the previous 
study' which showed that the calculations based on the 
one-dimensional heat flow agreed well with those based 
on two-dimensional heat flow. Hence the present 
design charts could be used to obtain the results of the 
previous study! in a much shorter period of time. 


Example 2 


Consider the wedge-shaped wing to be made from 
four different metals, such as aluminum alloy, beryl- 
lium, steel, and titanium. Initially the wing is flying 
steadily at a Mach Number of 1.4 at an altitude of 
50,000 ft. and then is accelerated at a rate of 2 Mach 
Numbers per min. to a Mach Number of 6 and finally 
flown steadily at 6. At the point in question along 
the chord, the half-thickness is 0.15 ft., and the local 
heat-transfer coefficient is given approximately by 


h = 36 (M — 1.4) 


Calculate the local temperature distribution when the 
wing just reaches a Mach Number of 6. 

Table 2 summarizes the properties of the different 
metals used and also the values of C; of Eq. (40). Note 
that Eq. (40) reduces to a single term in the final solu- 
tion for this particular example, since Cy is zero. This 
simplification is a result of the particular linear form 
used above for the approximate relation between /: and 
M. 

Table 3 summarizes the calculation of the tempera- 
ture distributions for the four materials, using the 
scheme described in detail in Example 1. The com- 
parison of the temperature distributions is shown in 
Fig. 9. The important fact to be noticed in Fig. 9 is 
that the engineering results can be obtained quickly 
from the design charts at any instant of time without 
having to solve the entire history of the transient tem 


perature. 


CONCLUSIONS 


Design charts have been prepared to permit rapid 
and reasonable estimates of the transient temperature 
distribution for many practical situations arising from 
aerodynamic heating at supersonic speeds. These 
design charts were prepared from the analytical solu- 
tion for the one-dimensional heat flow in a solid body of 
constant properties when the local surface coefficient 
of heat transfer is a strong function of time and when 
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Gas Dvnamic Behavior of Real Gases’ 


L. N. TAOt 


Worthington Corporation 


SUMMARY 


Gas dynamic problems of real gases can be treated in a most 
general way by the use of the Beattie-Bridgeman equation of 
state and a specific heat depending upon temperature and pres 
sure based on the spectroscopic measurements. The results in 
the 
geometrical quantities determining the flow configuration may 


form of gas dynamic relations between gas properties and 


be presented in terms of specified constants with the Crocco 


Number as independent variable 


(I) INTRODUCTION 


I THE TREATMENT OF GAS dynamic problems it 1s 
often assumed that the behavior of the medium is 
that of a perfect gas. In order to improve the accu- 
racy of the solution for real gases, several investigators 
have treated the relations of the isentropic flow and the 
Rankine-Hugoniot flow through plane shock by means 
of various assumptions. By considering deviations 
from both definitions of a perfect gas, 
Eggers” have obtained the results for one-dimensional 
flow problems based on Van der Waals’ and Berthelot’s 
with the 


Tsien! and 


equations, respectively. Dealing variable 
speciic heat, Tsien suggested a quadratic function of 
temperature for air only, and Eggers used an equation 
suitable for diatomic While neglecting the 
deviation due to the thermal imperfectness, Durham* 
derived results for air with the Boyle-Charles equation 
of state and the mean value of the variable specific 
heat, and Korst* supplemented the result in dealing 


with the problems of two-dimensional shock configur- 


gases. 


ation. Other attempts were based upon approximat- 
ing special processes for a real gas: Walker?’ analyzed 


the problem on the discharge of gases through an ori- 
~ : . : 1 AT 
fice or nozzle employing an adiabatic equation pve 


constant instead of pv‘ = constant. Traupel* 
attacked the problems of a real gas by assuming that 
the compressibility factor of a gas can be considered as 
a constant for constant entropy, but this concept can- 
not be applied generally, since not all gases lend them- 
selves to such an approximation. 

In the present investigation for improving the ac- 
curacy in gas dynamic problems of the isentropic flow 
and normal and oblique shocks of real gases, a more 
generalized concept is introduced so that it may be 
applied to any gases. The Beattie-Bridgeman equa- 
tion of state? is used as the thermodynamic relation of 
pressure, temperature, and density. Values of the 
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specific heat as a function of temperature and pressure 
are obtained from spectroscopic data and the use of 
the Beattie-Bridgeman equation of state and repre- 
sented in a descending polynomial. The gas dynamic 
relations are in a generalized manner in terms of the 
specific heat equation; these then may be applied to 
any gas after introducing proper numerical values for 
that gas. 


(II) BEHAVIOR OF REAL GASES 


Based on classical thermodynamics the real gas 
differs from an ideal gas in two significant ways, known 
as the thermal imperfectness and caloric imperfectness. 


(A) The Equation of State 


It is recognized that the simple equation of state for 
a perfect gas is only a first approximation for the be 
havior of the actual gas. Hence many modifications 
have been suggested; more than 50 equations® have 
been proposed to represent the actual gaseous behavior 
with closer agreement. Up to now it seems impossible 
to express the complete behavior of a substance over 
the entire range of the measured pressure and tem 
perature with a simple equation. 

Due to the impossibility of establishing a complete 
equation through the whole ranges, the empirical or 
semi-empirical equations arose. In the engineering 
field the most widely used equation by researchers is 


the Beattie-Bridgeman‘ 


[(1/p) + B] — Ap? (1) 


pb = RTp*(1 — « 

with 
e = cp/T* 
B = B,(1 — bp 


and 
A = A,(l — ap 


This equation gives satisfactory results for the thermo 
dynamic properties of real gases and fluids. This ex- 
pression for the pressure-volume-temperature relation 
ship over a wide range of temperatures and pressures 
gives 0.1 per cent to 0.2 per cent deviation for most 
gases. Hence this equation is used extensively here as 
the equation of state. The constants in the Beattie 
Bridgeman equation have been determined for a num- 
ber of common gases. Even though this equation is 
used for some uncommon gases and their numerical 
constants are not available, the constants may be 
evaluated approximately from the critical data by the 


suggestion of Maron and Turnbull’ with the assumption 





764 JOURNAL OF THE AERONAUTICAL SCIENCES 


that the activity coefficient for most gases is the same 
at the same reduced pressure and temperature. 


(B) The Specific Heats 

To supplement the purely thermodynamic methods 
of finding the specific heats, in most cases superseded 
by methods based on the statistical molecular theory, 
the required molecular properties of gases may be 
derived from the analysis of their spectra. The spec- 
tral data give essentially the various possible energy 
levels of the given kind of molecule, and the statistical 
molecular theory provides a general law of the distri- 
bution of molecules among energy states. By the 
introduction of certain generalizations taken from quan- 
tum mechanics, the statistical method may be applied 
satisfactorily even when lacking detailed knowledge of 
the spectrum. These theoretical values of C, and C, 
are in excellent agreement with the best known ex- 
perimental data. 

As suggested by Kelly" it is possible to represent the 
molal heat capacities of most gases empirically over a 
wide range of temperature by means of three-constant 
equations of the form 


C,, = i + (a2/T) + (a3/T?) (2) 


Such an equation, with empirically determined con- 
stants, can be fitted to the experimental data in many 
cases with a precision of one per cent over a range from 
273°K. up to 1,000°K. Using the thermodynamic 
identities 


(0C,/0p)r = —T(0%”/O0T”), 


(0C,/dv)r = T(0*p/OT”), (3) 

and integrating at constant temperature, we find 
C, = C,, — ST (0%/dT?), dp (4) 
C, = Cy, + ST (0°%p/0T?), dv (5) 


where C,, and C,, are functions of temperature only, 
and physically they are the corresponding values at 
zero pressure. Their relationship may be expressed as 
C.-C, 28 

and the integral parts of Eq. (5) may be determined by 
introducing the equation of state. Here applying the 
Beattie-Bridgeman equation of state, it becomes 
Cy = Coo + [(2A0h)/(RT?)] + (12cp/T*) + 

[((Dp?) /(RT?*)] — [(8£6%)/(R°T*)] (6) 


C, = Co + flv, T) = Cw + (6cpR/T*) (7) 


(III) ANALYSIS 


(A) Isentropic Flow 
For the isentropic case from the basic equations of 
fluid flow we get 


CaT + T(dp/oT), d(1/p) = 0 


or CT — (T/p*) (0p/OT), dp = 0 (8) 
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The Beattie-Bridgeman equation is introduced in the 
form 
p = pRT [1 — (cp/T*)] (1 + Bop — Bobp?) — 
Ay p*( ] ap 
or 
= pRT[1 — (cp/T*)] (1 + Bop — Dp?) — 

(Aop? — Ep*) (9) 
with D = Bob and E = Aya. Hence (0p/O07T), may 
be derived from this equation and substituted into 
Eq. (8) 

C.dT — (RT/p) (1 + Bop — Dp?) > 
[1 + (2cp/T*)]dp = 0 (10 


Also, as discussed in the previous section we may write 


C= C,, + f(v, T) = a, + (ae 'T) + 
(a3/T*) + (6cpR/T*) 


II 


Cif{1 + a + (B/T) + (y/T?) + 
(6cpR/C,T*)] (11) 


where C; is a constant that has the same value as the 
reference constant specific heat. With this, Eq. (11) 
becomes 


Cyl + a + (B/T) + (y/T?) + 
[(6cpR)/(C;T*) |{dT — (RT/p) X 
(1 + Bop — Dp’) {1 + (2cp T*) |dp => 0 (12) 


In classical thermodynamics any thermodynamic prop- 
erty of a pure substance may be expressed explicitly as 
a function of two other properties—e.g., T = T(p, S). 
In the present case where the entropy is constant, it 
follows that for a specified entropy any property may 
be expressed explicitly by another single variable, 7 = 
T(p). 

Since all the terms inside the braces and brackets of 
Eq. (12) are rather small quantities in comparison with 
unity, the well known “‘perturbation approach”’ method 
is now applied. It is assumed that all the high-order 
terms, which actually appear in dimensionless quan- 
tities, are neglected, so that the solution of Eq. (12) 
presents itself as 


T(p) = To(p) + aT ,(p) + BT3(p) + 
yTy(p) + BoTn,(p) + DT p(p) + eT Ap) (13 
This is substituted into Eq. (12). Since all the co- 
efficients are independent of each other, the following 
set of equations must hold true. 
CdT, — (R/p)Todp = 0 (14 
Ci(dT, + dT,) — (R/p)T,dp = 0 (15a 
Ci{dTg + (dT/To)] — (R/p)Tsip = 0 (15b 
C,'dT, + (dT); To?)] — (R/p)T,dp = 0 (le 
CdTs, — (R/p) (Ts, + pTo)dp = 0 (15d 
CilTp — (R/p) (Tp — p?T)dp = 0 (15e 
C.dT. + (6pR/T,*)dT — 
(R/p) [T. + (2p/Ty?) ]dp = 0 (158) 


These linear equations are solved and Eq. (13) becomes 
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T = \p'-! — ad(k — 1)p*—' log p + B+ (¥/(2Ap*—')] + Bok — 1)Ap* — 
D[{(k — 1)/2]Ap*t+! + C(2/d?) (Rk — 1) Dp®-* $= (16 


where \ is an integration constant and k = 1 + (R/C;). With the initial condition 7), and p) the complete solution 1s 


7 ) k—1 8 1—k 2—2k 
ype ome 2h-C"}-ab 0)" 1. 
;; p { pi TL po 2T, Po 
k— 1 2 I pp 13k] 
Byp(k — D(1 -") + Dp? [ -(*)'|- 2 [ (* \ (17 
Po 2 Po T,? p f 
To represent p, py in terms of 7° Ty, we get 
0 | ie a ca l r 8 l i ie l i 
LY peat te tsL CY] + haat ilt CY T+ 
ae t~-t "Rh" Re= 1 i i” Hea 1 Ty 
TY Ne—t ] D pv? (7) h | =| (=) 4$—3k)/(A j 
Bopo | 1 — | — | — 1 — (18) 
. (7) J 2 Ty " Ty T) iT 
Also with the aid of the equation of state, we get the following equations 
Yh vvet 26 O Rl-C 
. 41 — afk — 1) og — — b+ 7 j= “ 
ie os | p Ty? 2 po 
p ,kt+ 1 p \* po. teal 
Bopok | 1 — + Dpy’ r= 1 — ae 
po 2 po 3 " 
o— Be \ 3-A 
0-0") « 
Po RT, Po f 


p p\*( | k— 1 p B 1 | (7) . | yy 7 | (2) _ 
. 41+ log en ea a est 4 
p () was ke” py bi by "Te4 pi 
p " ad (2)" “po: 
Bopo | 1 — — Dp’ 1 — 
| (5) ‘ 2k Po ) Pr 
Apo | | (? _— 1 Ep? | | (*.) of ‘| 
ara x T° ae = 20 
RT, k Po RT, k \ po 


Similarly, 


T (*)" nik a a p 8 | (err) “ 4 (nen) 
—_ = <l—a : log —. 1]—{— — a = 
T, Po \ k? Po To k Po T,? 2k h, 
Dp? k — 1 [ 7 (2)") Coy 3k — 1 [ (2) 1—k " 
2 k Po T,? k Po 
k — 1 Aopo pn ] Ep? k — 1 | (*) 3—k ‘| 
= Tie > (22 
k RT, [ ( TRI, po ; = 


In these equations if all the correction coefficients that exist in dimensionless quantity forms are zero, then those 
equations are reduced to those for a perfect gas. 

For the explanation of such an isentropic flow, the characteristics and properties of a ‘Laval nozzle”’ are analyzed. 
First, from the energy principle and Eq. (18) we obtain 


udu + d(p/p) + {(pCdT)/[T(ap/oT),]} = 0 (23) 


~ 


With the expressions for (0p/07T), and C,, this equation becomes 


udu + d(p/p) + Ci} 1 + a+ (8/T) + (7 T?) — 3(3 — 2k) (Cp/T*) — [(Aop)/(RT)] + [(Ep?) (RT) |} dT = 0 
(24) 
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Integration of this equation with the initial condition vu = Oat 7 = 7) gives 
Catered) eo(-1) far fl) 
a= — TC ‘7 eS x c _— = aR log a _~- _ 
: . a” “ihe wy ti? T, 
3C 0) r (3—2k)/(k—1) Ay 7) ‘if Wc ; k _ l E 07 T 2/(h ; 
—- @-o[1-(G.) |-@-» in [1 - (7) |+ - [1 -(F L095 
T,? L 7 RT, T 2. R™ML T 


Now we introduce the so-called ‘‘Croceo Number," which is defined as the ratio of the local velocity to the ‘“‘Lava] 


velocity’''' or the maximum velocity attainable in the same direction. In the one-dimensional problems it is 
W = u/Umar = u/[2f$cp1T]"? = u/\u? + [2/(k — 1) Ju?) (26 
and its relation to the Mach Number is 
M? = [2/(k — 1)] [(W?/(1 — W?)] (27 
or W? = 1 — (1/}1 + [(k — 1)/2]4P}) (28 


Such a variable is used because it transforms the infinite range of Mach Number to a range between zero and unity. 
Also it is rather widely spaced at the low speed, and hence, it gives less calculation error. Furthermore, it is simpler 
to handle in the derivation of formulas. 

Before introducing the Crocco Number, the sonic velocity should be found. By the formula u, = (dp/dp)'/* we find 


, ap bo (p\*'s — ( ) | : (°) | 
tn wh i EY 33 eee - 
a dp Po (4) | ' k T = Po T k Po 

walt gE) fo aae(r e822) + aml 52 2) | 
an © | = _— By wR l D 0° l— 
Ty? 2 | k pi r ba k po "o k 4 | 

> po } —_ 2k i—38 41 Po 2 2—h % n- , 7 

3 — 2) [1 - (° ) |+aefi- (*) |-% |! - (* f (29 

T,? k Po RT, | k \p RT) k \p f 


The local sonic velocity is obtained simply by moving the reference point to the local state and is evaluated at this 
reference state, by) = p, p = p,and 7) = T. Soit becomes 


aL k— 1 Bk-1 7 i=}, ,k+ 1 cp (3 — 2k) (4 — 3k) 
== ’ 


2 =k! _ _ + Bop — D 
; p ke T k Tf k Co k 


Aop ( _ ‘) Ep? (’ -k 
aa - _— (30 
RT\ k RT\ k 
With Eqs. (24), (25), (26), and (30), a relation between temperature and Crocco Number is found. 


is T = (1 —_ W?) ‘1 + aah; + (B/T))te + (Y 'T 0?) ts + Bo pots + D po "ts a 
(cpo/To*)ts + [(Aopo)/ (RT) tz + [CE po?) /(RT>) ]ts} (31) 





For abbreviations let us put these dimensionless ex- tr = [(1 — k)/k] + (2/k) A — W2)2@-*/E-» 4 
pressions a = 2%, B/Ty) = 2, ¥/To” = 23, Bopo = 2, fl — (3/k)] A -— wee 
Dopo? = 25, Cpo/10* = 26, (Aopo)/(RTo) = 27, and (Epo?) + tg = [(k — 1)/2k] — (3/k) 1 — W2) 8-1-4 
(RT) = 23. Then [7/2k] — (1/2)} (a — W224 
r - ne = " By the substitution of Eq. (31) into Eqs. (21) and (18 
a? (: T p> sti) (32) the following equations are obtained. 
where p = (1 — wees ( 1 ¥ sin) a 
h = W? " = 
lp = [(k — 1)/k] (W2/1 — W?] - where 
(1/k) loz (1 — TV?) wm = [k/(k — 1)]t, + [1/(Rk — 1)] log (1 — W? 
ts = (1/k) (Wk — W2))/ — Wy? mw, = [k/(k — 1)]t — [1/(k — 1)] [W?/1 — W?)] 
ts = [(R — 1)/k] — 201 — W?) V@—-D 4+ m3 = [k/(Rk — 1)]fs + [1/(R — 1)] X 
[((k + 1)/k] 1 — W2)te@- {1 — [1/1 — W?)*}} 
fj = [1 — k)/k] + [2+ C1/R)] I — W)2/e-) wy = [R/(R — 1) ]ly 
(1 + (2/k)] 1 — W2)@+vd/e-1 mw, = [k/(k — 1)]ts + (1/2) — (1/2) (A — W?)*! 
te = }[4(1 — k)]/R} + ([(3 — 2k) (4 — 3k)]/R} X we = [k/(k — 1)]t¢ +3 — 311 — W2)4-8H/0-1 
(1 — W2)4-39/4—D — }1(4 — 3k) (2 — R)]/k} &X wi = [(R/(R-—1)]Jp +1—- (0 — W2)2-%/e-1 


(1 — W?)G-*)1e-1 TT, = [k (kh — lye —1 +1 — W2)8-¥/4—-D 








and 


wher 


As th 
cal Vi 
gas fi 
the c 
be ret 

Fur 
flow t 
the v 
They 
tinuit 


For ¢c 
taken 
by al 
pe int 


At thi 
ergy ¢ 
By st 
Eg. (3 
A 

At 


where 


Se 





aval 


(26 


9 


his 





GAS DYNAMIC BEHAVIOR OF REAL GASES 767 


and 


P= (1 — Wye "(1+ 2 sai] (34) 
1 l 


d) [4 /(k — 1) ] + [1/(k - 1)] log (1 — W?) 
[to/(k — 1)] — [1/(k — 1)] [W270 — W?)] 
d; = [ts/(R — 1)] + }1/[2(k — 1) X 


}1 — 1 — W?)?}} 
- ly (kR —1l)J+1—-Q - WwW? 
ds (t5/(R — 1)] — (1/2) + (1/2) 1 — W?2)2/¢ 
di [t6/(R — 1)] + 2 — 201 — W?)G-%)/-1 


d; = t; (k — 1) 
ds ts (Rk — 1) 


As the correction factors are evaluated with the numeri- 
cal values for a specific gas, the deviation of the real 
gas from the perfect gas can be determined. If all 
the correction coefficients vanish, these equations will 
be reduced to the well-known formulas for perfect gas. 

Furthermore, in the one-dimensional isentropic 
flow through a Laval nozzle the cross-sectional area and 
the wall force function also play important roles. 
They are analyzed in the following: From the con- 
tinuity equation we get 


d pu = (0 (3.9) 


For convenience in application the reference point is 
taken from the critical cross section, which is indicated 
by an asterisk instead of the previous stagnation 
point. Thus 

puA = p*u*A* (36a) 


A/A* = (p*/po) (po/p) (W*/W) (36b) 


At this particular section it is easily verified from the en- 
ergy equation that IV’* is equal to [(k — 1), (k + 1)]'”*. 
By substitutions of p*/p, and p/p into Eq. (36b), 


Eq. (36b) is solved. 
A 1 E= 19 | "(G3)" x 
= ( ~ - 
A* wl 2 k+1 
( + Y «S:) (37) 
1 l 


S, = [1/(k + 1)] + [1/(k — 1)] log [2/(k + 2)] - 
[t:/(k — 1)] — [1/(Rk — 1)] log (1 — W?) 


where 


S. = —[1/(2k)] — }1/[k(k — 1)]f X 
log [2/(k + 2)] — [t/(k — 1)] + [1/(k — 1)] 
(Ww2/ — W?)] 


Ss; = [(k2 + 1)/(4k)] — (Rk + 1)2/[8(R — 1)]} — 
[ts (k—1)]+)01 — W?) 2/[2(k — 1)]}} 
S; = [(k + 1)/(2k)] + [2/(e + 1)]'*-' X 
[— (k + 2)/k] — [4/(R — 1] —( — W276 
Ss; = —(1/k) + [2/(k + 1)]2/4-) X 
1[(k + 2) (k + 3)]/[2k(R + 1)]f - 
(ts/(R — 1)] — (1/2) (1 — W2)2/4-» 


Se = —(4/k) + [2/(k + 1)] *-) (3 — 2k) X 
()(4 — 3k)/[R(R + 1)) + [1/(R + 1)]}}) - 
[t,/(Rk — 1)] +201 — W230 


S; = —(1/k) + [2 (eo Jc? KIKI) ¥¢ 
(4/[k(R + 1)]} — [nk -— 1] 
Se = [1 (2k) | = [2 (R + 1)]° ky /(k—-1 x 


f4/[k(k + DI} — [ek — 1)] 


The wall force function F is defined as 


os 


F = pA + pu°A 
Hence, expressed in dimensionless form it becomes 


F/F* = [(p/po)/(p*/po)] (A/A*) 
(11 + (pu2/p)]/{1 + (p*u*?/p*)]} (39) 


Then with the aid of the equation of state we get 


pu? 2kW? [ bm 
p (k-1)(1— W) "OR 


8 (' -- , ‘ - ) " 
T\ k rm\ ere] 


»k+1 cp (3 — 2k) (4 — 3k 
k T? k 


Ao 2 Ep? [3 
”( -1)4 “(-1)| (40) 
RT \k RT \k 


Therefore Eq. (39) becomes 


Dp 


F (k+1)W2?+ (k-1 ( 2 ) 
— = : 1+ > 3.j (41) 
Fr 2W(k? — 1)!” i=1 
where 
1s I ' 2 
= = og — 
vai ' k—-1 -kR+1 k+1 
2(k l a! 


2 
es ee log - 
pines ae OS ee FS 

0k — 1)W? 


(k+1)W? + (kR — 1) (1 — W?) 


k—1 
j33 = 73 S3 + vin 
~N 
2(k — 1)W? | 
(k + 1)W? — (kR — 1) (1 — W?)? 
e+ 3/ 2 t 
= S; -— 1+ + 
)4 ™4 + { +1 are 
kW? - 
(1 — Wve 


(k + 1)IV? + (k — 1) 


ot FORT 2. ee 
| ll ae te 7 eT ke+ 1] 7 


2(k + 1)W? 


(k+ 1)W? + (k 


| 
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jo = te + So t+1 + where subscripts | and 2 refer to the upstream and 
[ 2(4 — 3k) (3 — 2k) W? ; downstream flows, respectively. By introducing the 
3— — 72) (4—3kh 1 . . : 
ig (k + 1)W? + (k — 1) (1 W*) Crocco Number, they become 


II 
4 


+5 -( | )( ; y" i ‘= pW, — pili = 0 (42 
iis k+i/\e +1 


2(2 — k)W? Ee pz — Pi + Um*(p2IV2? — pli?) = O (43b 
(1 — W?)@-*)/a-1 


(k + 1)W? + (k — 1) . BY Cp 
folrtati+ Zt oe - ) ear + 
' LS 4 | . 2(3 — k)W? 1 is lig 
Ig = We T Os T rs 
= Z (k -+- l VW % + (k = 1) pe Pr 4 3cR p, 3cR po { 
ire = 8 = RR PI) 
(1 — W2)6-H1e-1 p» p1 T;? T.? pP Pl 
(B) The Normal Shock E | y 
ar (po? — pi?) + . Un (We? — W,?) = 0 | 
The continuity equation, the momentum equation, the ” 2 ; es 
and the energy equation for the one-dimensional shock 
are The ‘‘Laval”’ velocity can be evaluated from either one 
polls — pity = O (42a) of the following equations 
pos” —_ piluy- + Ny - 5 = 0 (43a) Un? : » SC, IT 
ae ol I 
fcar+ | a(-) + ot 
. . Ov /7 p 
p2 Pr l ‘ ( 2 ) 
— + - (%* — 4,7) = 0 (44a) = a," > (45 
2 k—-11-W: | 


p2 Pl - 


Here the second form is used, since it can be expressed in terms of the properties and the Crocco Number 


ttm? 


_ Qk | a k k—1 _e=1 k+1 
—k-11-Wepl | 7 ; 


Then the three equations with the substitution of ‘‘Laval’’ velocity and the equation of state have four unknowns 
p», T2, po, and W.. Although it is a very tedious procedure, by the elimination of IV2, p2, and 72, they yield 


k 7 k pr (: = a) . 1— W,? p2 (1 7 =) [1 — (po/pi)]? (Rk — 1) (1 — W,?) 
k — | k — | Pi Pi 2W,? Pi Pi tkW,? 


Pa (k — 1) (1 — Wy?) % ( 7 an (: 7 a ‘ (k — 1) (1 — W.2) ‘ad 7 (: a ai ‘i 
al 2kW 2 pi Pi k—1 pili 2kW 2 Pi py 
4k?W,? f | 2kW,? pi pb} k-1 OVAL pi 
(kR-1) (1 —- m= (1 _ al (k-— 1)? (1 — W1*)\ eS (k — 1) 1 — W,2) at 7 P) 7 
2kW,? pr 1k2W2 f ™ % 2RW 2 pr pr 
{ 


— 
a, 2kW)? ( y’ (k-12(1-W2)) . 
sok |. + 2,(0) + 


= , l rs a 
P= 21 [1 — (bo/pi)] (k — 1) 1 — Wi) pi kV? 


xX 


[AB (1 -} RW jae 
~LkR—1 Pp (ek — 1) (1 — W,2) [1 — (po/pr)] + 2kW2S 22 py Pi 
J 2k" iN _b [1 — (p2/p:)] (k — 1) (1 — W2) 
\(k — 1) (1 _ W,’) [1 ~— (pe pr) ] + 2kW;? | Pi (k — 1) (1 a W,?) {1 —_ (po pi)] ne 2kwW 2 
[1 — (po/pi) |? (R-— 1) A a] | k pe p1\? 
4k? ide k—1p, (: a. (7) " 


f. 2kWy2 4 ) 1- WwW? ( 6 be — 28) (4-38) 
(ek — 1) (1 — WY) [1 — (p2/p)] + 2kW2h OW? pi \ k a 


a kh \" (*)\ + (1 2 *) _ fl = Gyp) PA - DG - Wr) 
(ek — 1) (1 — W,2) [1 — (po/pr)] + 2kW 1? ps pi 1kW 2 
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[(3 — 2k) (4 — 3k) | j 2kW,? (*\) ( (*) 
—1}/—41- > ¥ “eee x 
L k ' (k — 1) (1 — W,?) [1 — (pe bi)] + 2kW,? \po f P2 
2kW;? ) he 2kW? | : 
\(k — 1) (1 — W,2) [1 — (po’pr)] + 2kW 2S | (k — 1) (1 — W,?) [1 — (po/pi)] + 2kW25 
- [ k Do (1 Pis 2kW,? }) l — W,? (: alt o 
"tk— lb p(k — 1) (1 — W,2) [1 — (po/pi)] + 2kW2S OW? pJ\k 7 


bp. \(k — 1) (1 — W,2) [1 — (po/pi)] + 2kW 2h Pi 


2(k — 1) (1 — W,?) [1 — (b2/h))] [1 — (po'p,) |? (R — 197 1 - a) | k pr, 
(k — 1) (1 — W,2) [1 — (po/pi)] + 2kW 2 2k? W 2 ~ LI—k py 


J, py, 2kW;,? ( l - W, (: ) 3 6S 
ps (k — 1) (1 — W2) [1 — (po/p)] + 2kWES We a 


) 2kW;,? (2) = Poe ’) P 
\(k — 1) (1 — Wi?) [1 — (po/pr)] + 2kWi2) \pe pr 
3 2kW o\ (Rk — 1) (38 — 2k) (A — W,?) 
(3 ahi 2 4(1 A) G=DO—WO= MET 5 
2 \ l(k —-1l)(1 - W,?) [1 _ (p2 pi) | + kW, 24 Pi 4k?W," 


It is possible to simplify the solution of this equation by ignoring all the high-order terms of these dimensionless cor 


po (*) - (**) 
= + VA 
Pi pido p> Pils i - 


Since all the z,'s are independent of each other, each of them must equal 


rection factors, so that 


Substitute this solution into Eq. (47). 
zero. Equation (47) is now reduced to 


eta)" OL CO E-O)ta ab O-* 
b b= 1 Veh oW 2 pi/ Pil o Pil tk W,? Dio 


Its solution 1s 
( po Py , =] (50a 
or = [(k + 1)71V,? — (k 1)?] [(k? — 1) (1 — W,?)] (50b 


Ignoring the trivial solution, the nontrivial solution is the same as that derived for perfect gas. 
While these expressions involving Z; are solved and their solutions are substituted back into Eq. (48 


ps _ e+ tc al the “(i re, <4.) (51) 
Pi (k? — 1) (1 — W,?) t=1 
where 
k—-1 (kR-1)-( +1)? 
An's 
e+ 1 (e+ 1°We — (k — 1)? 


P (2-1/0 -WYwW? | 2 :; (k+ 13° (k — 1 | 
,_= a — og 
(k + 1)2W2 — (k — 1)? (R+ II -WY) (kR+DWP-(kR-1) ° (+120 - WYA)W? 


7 (2? -—1)(1— W2W? 5 { 2k + 12+ (k- 14) 
 (k+12We—-(kR-120 &+VYDAC-WY (+122 - (ek — 12S 
Ay = 0 
k—1)(1 — W,2) k+1 
A; = e — iE — = W,? + 
(k— 1) — (R + 1)W?2 UL k-—1! 
ow? |" + 1)2W,2 — (k — 1)? Ww? ) 
k(k — 1) (1 — W,?) (k? — 1) (1 — W,?) k(k — 1) (1 — W),21 
—_ (k? — 1) (1 — W,?) 2(4 — 3k)W?? J 2(3 — 2k) (k +1) (1 — W,?) 
a (k + 1)2W,2 — (k — 1)2 \ #—-Dad—-W, | (&+12W2—- & — ID? 


k+i1\2?__ +1... (k++ 1)Wi2+ (k -—1)) 
Ww W2+1] — 4k — 
(G5) Peer | (k + 1)2W,2 — (k — 1)2f 
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R-1-W) f 4-202 
A; : 
k+ 1)71V, (R—1)?{(R-1)a- WY, 
ke k+DU-W, Wy 
as e+ 12We—(k-1 (k+l(1— WW, 


This equation leads to the solution of the density ratio 


between upstream and downstream flow 


; k+l] 4 
. W,? ( ae > =i) (52 
Pp) k— 1] i=1 
where 
g G{(k + 1) — (k — 13?) /[2k(R — 1]; A 
+ \((k+ 1)Wi2 — (k — 1) /R} 
(lk + 1)? — (Rk — 1)7]/[2Rk(R — 
HA) + h((k + 1)Wi? -— (k — DIR 
g (((k + 1)? — (Rk — 1)7]/[2R(R -— 
1) i As) + (Rk + DW? — (Rk — 1) /R} 
(((k + 1)7142 — (k — 1)?]/[2R(R - 
1) Ay) — ,[(k 4 LW? — (k — 1))/(R = 1){ 
(k + 1)711," — (Rk 1)? 
£5 As + 
2k(k — | 
(k+1) (Rk + 1)W? -— (k — 1)] 
k(k — 1 
(k + 1)°W4" (k 1} 
ae 8 
2k(k — | 
(3 — 2k) (4 — 3k) [ (R+ 1)V2 - (k 2] 
} L k— 1 
(k + 1)71V,°7 — (Rk -— 1) : 
_ kk — 1 ~ 


3—-k[(k+D)W2-(k- | 
k b= 1 


With the pressure and density ratio equations the 
temperature ratio is then found with the help of the 
equation of state. It is 


ie | 


= [ (-—) a 
, t= k+1 W,2| 


where 
/ a 2 
h he ~% 
hs = As — g 
hy = Ay — gg +1 — [(R + 1)/(R — 1I)IWY? 
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directly from the continuity equation, 
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(C) The Oblique Shock 


In addition to the one-dimensional shock the oblique 
shock is analyzed here. Generally in such problems the 
velocity is resolved into two components—.e., the 
normal and tangential components with respect to the 
shock wave. Based on the previously stated concepts, 
the tangential velocity components do not change be- 
fore and after the shock, and the normal velocity com- 
ponents behave as in normal shock. 

Before we analyze any further the oblique shock, the 
Crocco Number in such flow should be discussed. Due 
to the addition of a velocity parallel to the shock wave 
front from the view of normal shock, the maximum 
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In addition to the one-dimensional shock the oblique 
shock is analyzed here. Generally in such problems the 
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normal and tangential components with respect to the 
shock wave. Based on the previously stated concepts, 
the tangential velocity components do not change be- 
fore and after the shock, and the normal velocity com- 
ponents behave as in normal shock. 

Before we analyze any further the oblique shock, the 
Crocco Number in such flow should be discussed. Due 
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tan (o — @), 1s defined for perfect gas and can be easily 
obtained from gas tables or a shock polar diagram. 
With the initial state and geometrical configuration, 
the true upstream Crocco Number may be found by 
this equation. Then the normal component of this 
number is calculated; so all the properties and char- 
acteristics of the oblique shock may be determined by 
the use of the normal shock equations in the preceding 


section. 


(IV) Discuss1on AND CONCLUSION 


As for numerical calculations of real gas behavior, 


correction factors of real gases are dependent upon the 
initial conditions. The formulas derived in the pre 
vious section have been applied to Freon-12 with the 
constants in the Beattie-Bridgeman equation, Eq. (1), 
determined by Buffington-Gilkey,'? as follows. 


Ae = 6.12 
By = 0.078 
a = 0.0404 
b = 0.0824 
c=0 


Values of P are in Ibs. per sq.in.; values of p are in Ibs. 
per cu.ft.; and values of 7 in °R. The constants 
in the specific heat equation, Eq. (9), determined from 
Raman _ spectroscopic measurements by _ Bradley” 
and infra-red spectroscopic measurements recently by 
Plyler-Benedict'* are as follows. 

a, = 26.768 

ado = —S8.3862 * 103 

a3 = 3.5825 X 105 


A series of calculations have been performed and tabu- 
lated with the Crocco Number as the independent 
variable.'® 

In order to compare the ideal and real gas behaviors, 
the following assumptions are made: The fictitious 
ideal gas has k = 1.14 and C, = 14.21 B.t.u. per Ib.- 
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mole-F, which are the numerical values at room con- 
ditions, while for the real gas the following initial con- 
ditions are chosen. 

(1) For isentropic flow the reference state was speci- 
fied to be the stagnation state, and the initial states 
were selected as (a) p = 30 Ibs. per sq.in. and 7 = 
500°R., (b) p = 10 atmospheres and 7 = 600°R., 
where the reference state of case (a) for the fictitious 
ideal gas applies with reasonable accuracy also for case 
(b). 

(2) In the calculations for shock relations, the state 
in front of the stationary shock was selected as the 
reference state and the initial states were chosen as 
(a) p = | atmosphere and 7 = 500°R. (b) p = 1/100 
atmosphere and 7° = 500°R. 

In the oblique shock problems, a wedge angle @ = 30 
was selected as an additional reference condition to the 
above listed initial states in both cases. In these ex- 
amples the initial states of lower pressure refer to the 
case that the real gas approaches the behavior of a 
semiperfect gas with only the caloric imperfectness 
being of appreciable influence, while the other ones are 
leading to results influenced by both imperfectnesses. 
All these calculations are presented in graphical forms, 
Figs. 2-11. 

In isentropic flow the pressure change deviates more 
from the ideal case than either the density or temper- 
ature change at the same state or Crocco Number, 
while in the shock problems, on the contrary, the 
pressure has less deviation. Also it is observed that 
the influences due to two imperfectnesses are generally 
of opposite direction so that they compensate each 
other partly. No numerical percentage of deviations 
are given, since these values will be different for dif- 
ferent assigned k and C, values of the ideal case. 

In this investigation additional phenomena should 
be discussed: First there is the heat capacity lag or the 
relaxation time required for the specific heat to change 
with a change of temperature. This time lag is mainly 
due to the fact that vibrational energy requires more 
collisions per molecule to attain a new temperature 
than the translational and rotational energies, and 
hence it reaches the equilibrium state after a longer 
interval. Kantrowitz'® found that if such time intervals 
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for the gas to absorb its full heat capacity due to the 


are comparable with or shorter than the time required 


change of temperature, the gas will not have reached 
its equilibrium state of energy. In this case the trans 
fer of energy from parts of the heat capacity which 
have more than their share to parts that have less 
than their share will be an irreversible process and 
will increase the entropy of the gas and consequently 
will decrease the total head of the gas. However, 
this may be made insignificant by lowering the rate of 
change of the temperature. As in the shock problems, 
Bethe and Teller'’ observed that the heat capacity 
lag will increase the thickness of the shock. This 
simply means a broadening of the region of transition 
As long as both the downstream and upstream positions 
are at a distance from the shock, this consideration 
may be ignored without great bearing on the accuracy 

Another possible effect that should be discussed is 
the chemical dissociation of the gas at high temper- 
ature. Such dissociation will lead to changes in the 
energy level of the gas; it will falsify the basic assump 
tion of the energy equation. But this can be avoided 
by suitable adjustment not allowing the gas to exceed 
its dissociation temperature. 

From the foregoing investigations the following con- 
clusions may be drawn. Gas dynamic problems can 
be treated in a most general way by the use of the 
Beattie-Bridgeman equation of state and a_ specific 
heat depending upon temperature and pressure based 
on the spectroscopic measurements. The results in 
the form of gas dynamic relations between the proper 
ties of the gas and geometrical quantities determining 
the flow configurations may be presented in terms of 
specified constants with the Crocco Number as inde- 
pendent variable. 
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T THE BEGINNING of the century, it was customary 
A to regard sound and fluid dynamics almost as two 
aspects of the same subject. The great names in this 
subject such as Rayleigh and Lamb were equally con- 
versant with both. With the rise of technical acoustics 
and aeronautics as separate subjects in which the ex- 
ponents of one group had little interest in the doings of 
the other, a dichotomy has taken place. The cleavage 
has been increased by the fact that the modern worker 
in acoustics requires an expert knowledge of electronics, 
whereas the instrumentation of aerodynamics has re- 
mained on the whole nonelectrical. Recent events, 
particularly the approach of aircraft speeds to that of 
sound and the increased speeds of ships, have, however, 
brought a new rapprochement between the two sub- 


jects. 


TURBULENCE AS AN ACOUSTIC PHENOMENON 


It is with Rayleigh that the story starts, in the era 
of Osborne Reynolds when scientists were trying to 
find the conditions under which a laminar ‘‘boundary 
layer’ (to use the modern term) in the flow over the 
walls of a pipe would become unstable and so lead to 
turbulence. Rayleigh’s idea was to take a steady 
gradient of velocity at the wall and impose on it a simple 
harmonic disturbance of finite amplitude to see whether 
the gradient in progressing remained stable. His 
mathematics showed that certain gradients were al- 
ways unstable, some neutral. Later Prandtl, von 
Karman, and Noether severally sought solutions of the 
boundary-layer equations for given boundary condi- 
tions in which they overlaid the steady stream with 
periodic fluctuations whose amplitude and wave- 
length depended on distance from the wall. 

The experimental method which this line of thought 
suggests is to impose on a steady stream a simple har- 
monic modulation of small amplitude and moderate 
frequency and to examine the transformation of the 
flow by this superposed fluctuation to see how closely 
it approximates in properties to a naturally turbulent 
stream, particularly in respect to (1) rise and decay of 
the fluctuations, and (2) development of characteristic 
distribution of velocity across the stream. The de- 
vices which have been employed to produce the de- 
sired simple harmonic motion in the course of the 
fluid comprise (1) tangential oscillation of part of a 
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smooth wall over which the fluid passes, (2) passage of 


the fluid over a rigid corrugated solid boundary, and 
(3) passage over a flexible boundary. The latter may 
take the form of a pennant or the fluid may be forced 
out of a linear orifice to impinge on a wedge so that the 
jet is set in sinusoidal oscillation between fluid bound- 
aries as in the familiar method of producing edge tones 
(see below). 

Most illuminating from the point of view of the 
architecture of turbulence are the experiments of the 
first type (Richardson and Tyler') in which a pipe has a 
loose section interposed between the entrance cone and 
the main portion of the pipe which is connected at its 
other end to an aspirator. The loose section is mounted 
on the end of a long connecting rod to a crank with ad 
justable throw on the axle of an electric motor. Air is 
first aspirated through the pipe at a speed just below 
the critical value and the gradient of velocity explored 
at various sections of the pipe downstream from the 
mobile portion. This is next set in operation and the 
velocity traverses redetermined. Whereas the former 
profile was parabolic, it has now gone over into the 
steeper logarithmic type characteristic of turbulent 
flow, which can in fact be reproduced in the untouched 
flow in the pipe by making the entry more abrupt, or 
copied at a larger scale by increasing the speed of as 
piration above the critical. Finally the distribution 
of velocity amplitude with the piston alone in opera- 
tion—i.e., no through draught—is obtained. The 
latter distribution can be deduced theoretically on the 
assumption of laminar flow. For a displacement wave 
of periodic form a sin wf having pulsatance w the veloc- 
ity amplitude aw reaches a maximum at a distance 6 
from the wall of the pipe, such that 


V w/2y = 2.28 


v being the kinematic viscosity of the fluid (see Fig. 1, 
curves at left). 

The disturbances of all frequencies are not equally 
amplified and some may be damped, though this de- 
pends also on the Reynolds Number of the flow. A 
re-examination of the earlier results showed a variation 
of amplification coefficient with frequency exhibited in 
Fig. 2, the broken line giving the theoretical variation 
(with square root of frequency). Such variations have 
also been shown in recent measurements by Schu- 
bauer and Skramstad and others ** 2” who also built 
loud-speaker diaphragms into the wall of a low tur- 
bulence wind tunnel with similar intention. 

It may be objected that this artificial turbulence does 
not correspond with nature in that the pulsatance w is 
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Fic. 1. Amplitude of velocity fluctuation in pipe of radius R; 
y, distance from wall. Left, in sound waves; right, in turbulent 


flow, u root mean square axial fluctuations, J’ mean stream 
velocity. 


not monochromatic in real turbulence but is more like a 
band spectrum embracing a range of frequencies re- 
quiring a Fourier series in the solution. To the extent 
that such partial vibrations may in practice have differ- 
ent damping coefficients so that the spectrum changes 
en marchant, our idealized picture of turbulent flow re- 
quires modification; but it is certainly a fact that nat- 
ural turbulence shows velocity profiles and distributions 
of intensity of fluctuation precisely similar to those in 
our developed flow. Thus Fage and Townend? in their 
ultra-microscopic investigation of turbulent flow in 
pipes have found peaks in the fluctuations near the 
walls like those shown, for the sake of comparison, on 
the right of Fig. 1. 

Although it seems then well established that tur- 
bulence in wind tunnels has its origin in a quasi-periodic 
disturbance, there seems to be no practical relationship 
between the frequency of the naturally induced oscil- 
lation and any other parameter of the flow, such as the 
Reynolds Number of the flow or the roughness of the 
walls. Any simplified periodicity that remains after 
a sufficiently long passage is to be ascribed to the greater 
damping of the high frequencies, which are, however, 
constantly replenished by new disturbances engendered 
at the walls. 

More is known about the onset of oscillation in the 
wake of an obstacle, such as a cylinder in a steady 
stream. Even at very low Reynolds Numbers a slight 
divagation of the wake is sufficient to set the wake 
pendulating just as though it were a free jet of liquid 
from a nozzle, set in motion by a slight vibration com- 
municated to the nozzle. At a somewhat higher 
Reynolds Number, the stationary vortices which are set 
up behind the obstacle take part in the oscillation of the 


wake and give rise to the familiar procession of vortices 
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down the wake, which, especially if the obstacle is Cap- 
able of sympathetic vibration, lead to the aeolian tones 
well known in acoustics. 

This motion is maintained by the disturbing action of 
the detachment of vortices from one half of the bound- 
ary layer on the pressure gradient along the other half. 
There are, as Holtzmark ef a/* have recently remarked. 
deviations from ideal flow in the boundary layer, 
leading to peaks of velocity near to the cylinder wall 
corresponding to those already remarked in pipe flow. 

All these oscillations, if of sufficient frequency and 
amplitude, give rise to sound. Which leads one to the 
question: Is there indeed any difference between such 
velocity fluctuations and more conventional sound 
sources? One necessity obviously is that the velocity 
changes should induce local pressure changes, to be then 
propagated with the speed of sound through the fluid. 
An instance of velocity waves not giving rise to sound 
waves is to be found in the propagation of shear vibra- 
tions. If a dise is oscillating about its axis in a vis- 
cous fluid it does not compress the fluid which bathes it 
and so it produces no sound. 


CORRELATIONS IN VELOCITY AND SOUND FIELDS 


In a turbulent fluid flow there exists a correlation 
in the statistical sense between the motions at two 
points in a field, but there is no fixed relationship (in 
time) such as one expects in the field of a sound source. 

One could say that correlations in a field of flow are 
vectorial whereas those in a sound field are rather of 
scalar type. 

An attempt to find a correlation between pressure 
fluctuations in a wind tunnel, picked up on a hot-wire 
microphone with streamlined nose, and the noise of the 
fan, picked up on a moving coil microphone, has been 
made by Schuh and Winter’ at the Royal Aeronautical 
Establishment at Farnborough. They analyzed the 
frequency spectrum of both types of disturbance, but 
were unable to establish a sure correlation between them. 
Nevertheless, it is recognized in aeronautical research 
establishments that in building a wind tunnel of low 
turbulence, noise is a factor which must be considered. 
When all other disturbing factors are eliminated, it is 
now possible to reduce the intensity of turbulence 
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Fic. 2. Amplification (+a) of periodic disturbances, propa- 
gated along a pipe, 5-in. diameter at mean speed 3 ft. per sec., 
in relation to frequency. Full line, experimental; broken line 
theoretical. 
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ACOUSTICS IN RELATION TO AERODYNAMICS iad 


expressed as a ratio of the fluctuation to the mean 
wind) to 1 3 per cent, which probably represents the 
fluctuation associated with the sound output of the 


propeller and other unavoidable sources. 


ACOUSTICS OF JETS 


Interest in jet noise has revived interest in another 
old problem, that of the “‘edge tone.” An unassisted 
jet produces an impure sound but may by some form of 
resonance or interference be persuaded to produce a 
series of pure tones. One way of doing this is to com- 
municate harmonic vibration to the nozzle from which 
the jet issues (Fig. 3a). Another way is to let an air 
or water jet impinge on a sharp wedge. If the dis- 
tance from slit to edge is kept constant and air at low 
pressure admitted to the wind chest behind the slit, the 
pitch of the edge tone rises with efflux velocity (with 
some jumps of frequency at certain pressures), while if 
the wind pressure is kept constant and the wedge moved 
away from the slit, an opposite effect on the pitch occurs. 
After the alternate vortex street behind an obstacle in a 
stream had been observed and studied, a similar vortex 
system was invoked to explain what could be seen in an 
edge-tone jet, although it was recognized that the eddies 
do not emerge full-grown from the slit but gather 
strength as they are carried toward the edge (see Fig. 
3b). 

About twenty years ago, two rival theories of the 
edge tone were current. One regarded the phenomenon 
as purely acoustic—each time a vortex hit the edge a 
sound wave signalled back to the edge to start another 
one. The author thought that he had disposed of this 
theory by opening another jet alongside the principal 
one and showing that the postulated sound wave is un- 
able to set this subsidiary jet into pendulation; but 
recently Powell’ and Nyborg’ have revived this acoustic 
hypothesis. The second theory was a hydrodynamic 
one, based on the writer’s* observation that a secondary 
vortex system is set up in the boundary layer alongside 
the wedge, and that these have to fall into line with the 
main system to form two vortex streets, one each side 
of the wedge. Curle® has recently taken up and ex- 
tended this theory, showing that the hydrodynamic 
reaction on the fluid debouching from the slit results 
from the alternating transverse velocity at the tip 
At low 
velocities, when the edge tone results from a simple 


without the intervention of any sound signal. 


periodic pendulation, it ought surely to be possible to 
give a purely hydrodynamic explanation to the phe- 
nomenon. It may be noted, indeed, that pendulations 
similar to those of edge tones may be produced when 
underwater jets strike an opposing wedge at speeds far 
too low to engender sound waves in the water.’ 
Reverting to the unassisted jet (from a circular ori- 
fice), one may assume that the origin of -the hissing 
noise which accompanies the efflux of air at moderate 
speeds is to be found in shearing forces, both within the 
nozzle and outside where a “‘surface of discontinuity” 
as Helmholtz called it—exists between the fast moving 
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Fic. 3. Jets in air made visible by admixture of smoke 
a) periodic disturbance communicated to nozzle, (b) jet imping 
ing on a wedge 


and the more or less stationary air. Here boundary 
layers engendering turbulence are to be found, at first 
apparent as small eddies or high-frequency fluctuations; 
but farther down the jet these roll up into larger ones, 
involving fluctuations of lower frequency. Examina- 
tion of the spectrum of open-air turbulence with hot 
wires and of the associated sound spectrum with micro- 
phones show them to be closely parallel. At any rate, 
confirmation of the statement about the relative 
size of the jet eddies can be found when one traces 
back the sound picked up from the microphone to its 
origins, the high-frequency sources then being revealed 
as located near the nozzle, while the low-frequency 
sources are to be found farther downstream and carried 
with the jet. 

When the speed of efflux from the nozzle exceeds that 
of sound in the gas, a change in the disposition of the 
jet ensues. Schlieren photographs show that the jet 
alternately expands and contracts, while it is criss- 
crossed with a system of shock waves originating at the 
nozzle and reflected at the confines of the jet. The 
phenomenon is now acoustic rather than hydro- 
dynamic in mechanism. One can, for example, let the 
supersonic jet—from a circular nozzle—impinge on a 
coaxial metal ring to produce edge tones of a definite 
frequency, but only if one places the round edge near 
one of the zones of instability at which the shock waves 
are reflected. Presumably in this case there is reso- 
nance of the waves reflected from the ring with those set 
up at the nozzle, as suggested by Wachsmuth for the 
low-speed edge tones. The jet-edge system just 
described is, in fact, the basis of the Hartmann gener- 
ator in which, by constructing the apparatus on a minia- 
ture scale (nozzle about 0.5 mm. in diameter, jet-edge 
separation about 1 mm.) this investigator was able to 
construct a source of ultrasonic frequency. 

A jet of supersonic velocity is, in fact an efficient 
transducer of kinetic into acoustic power, some of it 
audible, some ultrasonic. It has been suggested by 


Richards and others” '* ” that the origin of this 
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excess noise (above that of a merely turbulent jet) 
arises from the interaction of the shock waves with the 
turbulent boundary layer or ‘‘surface of discontinuity.”’ 
On this idea the zones of instability aforesaid will be- 
come potent sources of sound. The sound field around 
the jet will then be rather like that of light in front of an 
optical grating—that is to say, there will be directions 
in which the energy from different sources will con- 
structively and destructively interfere—and there is 
some experimental evidence, from polar curves of super- 
sonic jets, which supports this. 

These conclusions are fairly well supported by ex- 
periments on model jets. Such ‘smoothing’ of the 
polar curves for intensity in comparison with theory as 
is observed can probably be ascribed to the scattering 
of the sound as it passes through the turbulent field. 


THERMO-ACOUSTIC PHENOMENA 


Another phenomenon of respectable antiquity which 
has been resuscitated in this gas-turbine age is the ‘‘sing- 
ing flame.’ Actually, this was the first problem 
given to the author when in 1922 he went to the late 
Prof. C. H. Lees of London and said that he wanted 
to start research in acoustics. Indeed, the writer be- 
lieved that, as a result of his efforts, the matter had 
been decently buried, tut two symposia held at the 
Massachusetts Institute of Technology in 1952 and 
1954 have dealt inter alia with these thermo-acoustic 
phenomena. 

The sounds generated in jet engines and rocket tubes 
by the combustion of fuel at some point therein hark 
back to the Rijke tube experiment. This scientist 
found that when a heated gauze was placed in the 
lowcr half of a vertical open tube, the latter sounded its 
natural frequency vigorously as long as the gauze re- 
mained hotter than the air column. The interest to 
physicists in these oscillations lies in the peculiar man- 
ner in which they are maintained. It seems that the 
heat is in most cases applied to the gas for a small frac- 
tion only of the period of vibration to give a compres- 
sional impulse to the fluid locally. As long, however, 
as this impulse occurs at the right phase, the vibration 
of the column is maintained. In effect, as Rayleigh 
showed, this means that the source of heat must be 
located in the lower half of the tube, but the obverse 
phenomenon, ascribed to Boscha, in which a column of 
hot air is caused to sing by the application of cold in the 
upper section, has been suggested as the cause of noise 
in certain cooling towers of power plants. 

The relevance of these thermo-acoustic oscillations 
to the burning of pyrotechnics and rockets containing 
solid propellants was illustrated at the 1952 Symposium 
by Maxwell!! and by Smith and Springer,!? all of whom 
found instability and changes in the rate of combustion 
when the rocket started to ‘“‘whistle.’’ The latter pair 
of authors, verifying an idea put forward by Grad,'* 
suggest that the erratic combustion of aircraft rockets 
arises from the production of sounds which are usually 
one or another of the proper frequencies of the pipe 
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formed by the hollowed-out center or external annulys 
of the rocket charge. 


SOUND SCATTERING AND TURBULENCE 


The mode of attenuation of sound in the open air has 
become another matter of urgent interest since the 
noise level of aircraft propulsion has risen so alarmingly, 
Generally speaking in a fluid, apart from a fall of in- 
tensity from the source due to the increase in surface 
area of the wave front, the sound amplitude falls off 
with the square of the frequency, so that the ultrasonic 
and high sonic frequencies in a jet noise are fairly 
rapidly attenuated with distance. This effect is mag- 
nified often by the humidity of the atmosphere. 

The early investigators on the propagation of sound 
in the atmosphere used either subjective observation or 
objective apparatus incapable of following rapid 
fluctuations. Consequently they were obliged to 
explain their results in terms of the continuing meteor- 
ological elements like mean wind and temperature, 
though they recognized that certain meteorological 
conditions led to an irregular signal when picked up at a 
distance from a steady source. It is only recently 
with the introduction of sensitive electrical devices 
that studies of the fluctuations in the signal received 
from a steady source of sound have been possible. 

Another important effect is due to scattering of the 
sound by the eddies which exist in the atmosphere 
whenever it is turbulent—-which, near the ground, is 
almost always. This turbulence is, to within a few 
inches of the ground, isotropic, that is to say, the eddies 
are of symmetrical form and of size from many yards 
down to an inch or two. As such they scatter sound 
rather in the fashion that the globules of water in a 
fog scatter light. The relevant parameter is the rela- 
tion between circumference of the scattering entity 
(eddy) and the wave length. The air velocity in and 
around an eddy is of course fluctuating so that the 
scattering of sound is inconstant, rather like the flicker- 
ing of the light from stars due to atmospheric variations 
of density. 

These departures from the classical laws of sound 
propagation are to be ascribed to the turbulent nature 
of the atmosphere. The normal velocity c of sound in 
still air or in steady streaming is locally replaced by 
c + Ac, consequently we should expect a scattering of 
the sound from the direct path on the basis of the mean 
correlation coefficient R of turbulent motions over the 
path from source to listening post. Attempts have 
been made to set up theories which regard the eddy 
length as playing the part of the regions of abnormal 
refractive index which appear in Rayleigh’s theory of 
light scattering, with its dependence on wave length A 
and on the volume of the scattering entity in the form 
of v or /’, |. being the mean eddy diameter. There 1s, 
however, this difference, that the atmosphere is a 
medium of variable acoustic refrangibility without 
variable density—in the absence of thermal currents. 

Lighthill'* has applied his theory of sound generated 
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ACOUSTICS IN RELATION TO AERODYNAMICS rk 


aerodynamically to this problem and has shown that 


the energy scattered per unit time from unit volume is 
(S7u? ‘?)-(1T /c?) 


] being the incident sound intensity, / the mean eddy 
diameter, and u* the mean square velocity fluctuation. 

It may be objected that such scattering does not 
attenuate sound but merely redistribute it. While 
this is true for a point source, whenever the source is 
directional as most jet sounds are, at any rate for the 
higher frequencies, scattering tends to smooth out the 
polar curve of intensity distribution and to upset phase 
relationships between sounds arriving at the listener by 
different paths, so that in fact the decibel loss per 
1,000 ft. is greater in such an atmosphere than it would 
be in still fluid, the loss increasing with the square of 
the velocity fluctuation as Lighthill's formula shows. 

The author?’ has carried out measurements of velocity 
fluctuations, by hot wires; temperature fluctuations, by 
thermocouples; sound intensity at a distance up to 
100 ft. from a pure tone source, by microphones; and 
verified a number of features of the principles above 
expounded. In particular he has shown that some 
correlation exists between the fluctuations of sound in- 
tensity over the test range and intermediate fluctua- 
tions in both wind and temperature. Fig. 4 shows, for 
example, simultaneous fluctuations of velocity (as 
variations of hot-wire resistance) and sound intensity 
picked up at a height of 6 ft. over a grassy field. Such 
records, though they do not show an exact correla- 
tion between turbulence and sound propagation, do 
illustrate the effect of turbulent conditions on the root 
mean square sound pressure fluctuation within the 
medium, particularly when eddies comparable in size 
with the wave length of the sound are concerned. 


ACOUSTIC RELAXATION IN AERODYNAMICS 


Finally, reference should be made to some acoustic 
aspects of the use of freons and other vapors to replace 


air in supersonic wind tunnels. The raison d’étre of the 
substitution is, of course, to attain a given Mach Num- 
ber in the flow for a smaller expenditure of power, since 
these vapors have low values of the velocity of sound. 
They have different values also for the ratio of specific 
heats y but the users are prepared to correct for this and 
seem to be able to do so with success at least for Mach 
Numbers of the order of unity as a recent report by 
von Doenhoff, Braslow, and Schwartzberg?® shows. 

Certain questions arise, however, as to the exact cor- 
respondence of fields of flow even when the Reynolds 
Number and Mach Number are correct. The principal 
questions are: 

1) Does the change in y significantly alter the flow 
field and, if it does, can a correction be made for this? 

2) Is the equation of state of the same form? 

3) Will the vapor condense due to the lowering of 
temperature in an expansion? 

With regard to the variation of flow field with y, 
von Doenhoff, Braslow, and Schwartzberg!* have shown 





Fic. 4+. Simultaneous records of velocity fluctuations in open 
air (above) and sound intensity (below a, sound of 4,000 
cycles per sec., 6, sound of 7,500 cycles per sec 


that, at a particular speed, variation of y can be com 
pensated by a change in the relative thickness of the 
model although a change in y from 1.4 in air to 1.13 in 
freon involves an almost negligibly small change in the 
thickness ‘chord ratio. The speed range over which 
freon is an adequate substitute for air is, however, 
limited and is in fact from low speeds up to about 1.2 
times the speed of sound. 

If the relaxation time differs from that of air the 
whole flow pattern may be altered and some of the 
energy stays in the vibrational mode. The work of 
Miyahara and the author on the velocity of sound in 
the freons is pertinent. The relaxation time in freon 
12, for example, is 1.5 X 10~* sec. whereas that in air 
at sea level is 10~° sec. This means that a shock taking 
place in freon at a local JJ = 2 (U’ = 300 m. per sec.) 
at a distance of 1/20 mm. in front of the nose would be 
completely relaxed in this distance, whereas the shock 
front in air would be scarcely affected in the corre- 
sponding passage. This one would think must affect 
the local pressure on the model." 

Chapman'* who has recently examined the pre- 
sumed behavior of models in substitute vapor tunnels 
in relation to full-scale flight at high altitudes points 
out that the important parameter in this connection 
is not the relaxation time itself but the product of it 
and the thermal capacity of the molecules in vibration. 
He shows that thermal capacity lag should be important 
in boundary layers at low Reynolds Numbers or high 
Mach Numbers. It may be important at high alti- 
tudes, so that to this extent, a tunnel filled with a 
relaxing vapor or mixtures of vapors may be superior to 
one filled with atmospheric air, especially for hyper 
sonic flow. 

Proximity of conditions in the working fluid to the 


2 Assum- 


critical state affect both questions (2) and (3). 
ing, say, van der Waals’ equation, it would be possible 
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to make some modifications of the thermodynamic 
expressions for the pressures and temperatures in 
terms of the constants a and b of the equation, but a 
general formula of correction would be difficult to ob- 
tain. It has been estimated that the temperature 
may rise in places to a value of y~! times the static 
temperature. This sets a limit to the temperature at 
which a tunnel containing a vapor can be operated 
to avoid condensation. Also, there is a disadvantage in 
power consumption if one works near the critical state, 
for then the velocity and attenuation of sound are rap- 
idly changing. 

This rules out carbon tetrafluoride (one which has 
been suggested), but we are making sound measure- 
ments on the hexafluorides of sulphur, selenium and 
tellurium which look promising. 

What precisely the high absorption coefficient in 
these vapors (as compared to air) portends in this appli- 
cation is difficult to say but it must, one would think, 
have an effect on the devolution of the shock wave 


front. 
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The Ultimate Strength of Multiweb Box 
Beams in Pure Bending 


ROBERT A. NEEDHAM" 


University of California, Los Angeles 


ABSTRACT 


A method for predicting the ultimate bending strength of multi 
web box beams having formed sheet metal webs is presented 
Equations are derived on the assumption that the ultimate 
moment is reached when local crippling occurs at the junctions of 
the webs and the compression skin 
pared with 49 beams fabricated from 75S-T6 sheet and tested by 
the NACA. The 


tremely good 


The proposed theory is com- 


agreement between tests and theory was ex- 





NOTATION 


flange width of shear web, in 


Bb total width of box, in 
( = edge support coefficient for compression crippling, di- 
mensionless 

E = Young's modulus in compression, ksi 
E = secant modulus in compression, ksi 
E = tangent modulus in compression, ksi 
F compression yield stress (0.002 offset), ksi 
F = ().7 secant yield stress, ksi 
F compression crippling stress of skin, ksi 
F compression crippling stress of flange of shear web, ksi 
} = equivalent box height, in 
h = effective box height, in 
/ allowable bending moment of interior web unit, in.Ib. 
V7 exterior web correction, in.Ib 
nt = allowable moment per inch of box width, in.Ib. /in 
’ = number of webs in box 

= allowable load in compression flange, lb 

= skin thickness, in 

web thickness, in 
- = von Karman effective width, in 
\ = exponent in Ramberg-Osgood equation, dimensionless 
inelastic compression modulus ratio, dimensionless 
t = equivalent flange }/t, dimensionless 

ko = 1083V nWE/F., 


INTRODUCTION 


oo TRENDS IN WING DESIGN are toward the use 
of heavy skin surfaces employing various schemes 
of stabilizing the compression cover. The more com- 
mon types under investigation include multiweb, post- 
stiffening, double skin, and deep stringer constructions. 
Of the various types of structural configurations for the 
thin wing, the multiweb box consisting of heavy plate 
surfaces with a series of moderately light sheet metal 
webs appears promising, both from the standpoint of 
simplicity of construction and structural efficiency. In 
its simplest form, the multiweb box is shown in Fig. la. 
For the purpose of analysis, the wing structural box is 
idealized as shown in Fig. 1b. The variable height, /, of 
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the actual box is replaced by an ‘“‘equivalent”’ height, /, 
resulting in a rectangular box of the same width, B, and 
number of webs as the actual wing. 

The ultimate strength of a multiweb box beam in 
pure bending will be reached when the crippling stress 
at the junction of the webs and skin on the compression 
surface is reached. Failure may or may not be pre- 
ceded by skin buckling, depending upon the proportions 
and material of the box. Crushing of the webs normal 
to the plane of the skin (general instability) can be 
eliminated as a failure mode in a well-proportioned box 
since there is adequate experimental evidence to indi 
cate that the crushing forces in the webs are small.! 


BENDING STRENGTH OF MULTIWEB Box BEAMS 


Consider the rectangular box beam of Fig. 2a which 
has n vertical webs riveted to the upper and lower skins 
In the analysis that follows, the following is assumed. 


(1) Spanwise rivet spacing is such that inter-rivet 
buckling is precluded. 

(2) The distance }; (Fig. 2b) is a minimum con- 
sistent with good fabricating practice. 

(3) The effective bending element contributed by a 
vertical web is an equal angle having leg widths ) and a 
thickness /,,. 

(4) The effective bending element contributed by 
the skin consists of two effective widths of skin, w, 
having a thickness f,. 

(5) All vertical webs are of the same thickness, the 
upper and lower skins are the same thickness, and the 


flange width, 0, is constant for all flanges 
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Fic. 1. Actual and idealized wing box 
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Fic. 2. Interior bending element. 


On the basis of the above assumptions, the effective 
bending section of an interior web-skin combination is 
shown in Fig. 2c. 

In a recent paper® the author presented a method of 
predicting the compression crippling stress of formed 
sheet metal structural shapes. The proposed method 
was based on the assumption that a formed structural 
shape consisting of a series of flat plate elements could 
be treated as a series of angle sections possessing various 
degrees of edge support parallel to the direction of 
loading. A general equation of the form 


Fee = (CoV Fey) /(b'/t)” (1) 


was developed from extensive tests of channel- and 
angle-sections of 24S-T plain (unclad) and 24S-T3 and 
75S-T6 clad sheet. Fig. 3 illustrates the significance of 
C, and b’. Eq. (1) can be placed in convenient nondi- 
mensional form by plotting F../V F.,E against (b’/t). 
Such a plot, applicable to the alloys mentioned above, 
is shown in Fig. 4. 

Well-proportioned box beams of the type under con- 
sideration will fail in bending when the local crippling 
stress at the web-skin junction is reached. In isolated 
cases, the material of the web will be different from the 
material of the skin, and stresses in excess of the pro- 
portional limit of one of the materials of the combina- 
tion will result in a stress of a different magnitude in the 
other. A close approximation of the local crippling 
failure of these combinations can be obtained by as- 
suming equal strains at the point of attachment of the 
web and the skin. Consider the skin-web combination 
of Fig. 2b. The crippling stress, F,,,,, of the assumed 
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web angle can be obtained from Eq. (1). Since both 
flanges of the angle are supported (one flange by the 
web and the other by the skin riveting) the edge sup 
port coefficient, C,, of 0.366 1s used (Fig. 3c). The 
crippling stress in the sheet (= F,,,) at the web crippling 
stress is obtained from the combined compression 
stress-strain diagrams of the two materials, as shown jn 
Fig. 5. A “weighted” crippling stress of the skin-web 
combination is obtained in the usual manner from 


' = |( F cos Vir b + ( F.. tO. 1) te b + tw, yi 


Only the more usual case of identical materials in the 
webs and the skin will be treated here. At failure, the 
total axial load in the compression flange is given by 


ry = F...[n(2Wels T 2bt,,) — 2w,t, 
= 2w,tsFee}n|(bty)/(Wels) FI) — 14 (B 


where 2w, is the effective width of the skin acting at the 
web flange stress and m is the number of webs in the box 

For purposes of analysis, the structural index, as sug- 
gested by Shanley,’ will be introduced here. For the 
“wide column’’—i.e., the entire compression flange of 
the box—-the structural index becomes 1//Bh,*. Note 
that 1//B represents the moment per inch of chord, m. 
If the resisting moment at failure is defined by Pj, 
Eq. (3) becomes 
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m/h.? = [(2w.t;)/(Bh.) F..} n{ (btw) (wt,) +1) 4+ 1} 
(4) 


The effective width of sheet may be obtained from the 
von Karman equation‘ modified to account for inelastic 
effects. This equation is 


w/t = 0.85V 9 V E/F.e (3) 


Values of n are obtainable from Stowell’s Equation® for 
long plates simply supported on four edges—viz., 


n = (E,/E)[((1/2) + (1/2) V (1/4) + (3/4)(E,/E;)] 


(6) 
Substituting Eqs. (1) and (5) in Eq. (4) and simplifying 


m kN FE 
h (Bh, /t.2)(b/ty) 


where k = 1.03NV 9 VE F.,. 


COMPARISON BETWEEN THEORY AND TEST 


A well-proportioned aluminum alloy box beam of the 
type under consideration will develop a stress at the 
junction of the webs and the skin close to the yield 
stress of the material, in which case, E,/E = 0.7 and 
E/E = 0.2. The corresponding value of » [Eq. (6) ] 
is 0.59. In any case, however, 7 can be obtained from 
Eq. (6) at the web flange crippling stress, F,, = 0.366 
VV FyE/(b te)'/*, and the stress-strain diagram of the 
material. 

Reference 1 contains the results of 53 pure bending 
tests of 75S-T6 multiweb box beams of the type for 
which Eq. (7) was developed. Of the 53 tests, 49 had 
four vertical webs; hence, in Eq. (7), m = 4. The 
average value of F,, for these tests was 72,200 lbs. per 
sq.in. Assuming a value of Young’s modulus of 10,500 
ksi, Eq. (7) becomes 
m 3,100,000 V n 


he? (Bh, /ts?)(b/tw) 

{, | 0.204(b/t,,) ° ) 

4 ; - +1 “a (S) 
V 0 (ts/tw)? 

The tests of reference | were divided into three groups. 

Within each group the following quantities were held 


constant: 
Specimen Av. skin 
numbers b/t, jt. thickness, ¢, 
1-14 17.3 3.70 0.188 
15-35 14.4 2.44 0.125 
40-53 11.3 1.59 0.083 


To facilitate the determination of V7, Fig. 6 was 


prepared.* Since reference | did not contain compres- 


* See Appendix A for development of Fig. 6. 
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sive stress-strain diagrams, the Ramberg-Osgood equa- 
tion,® 

Fe/For = (F/Fo7z) + (3/7)(F/Fy.7)* 
with a value of x of 10, was assumed to hold. The 0,7 
secant yield stress was obtained from the equation’ 


For = (F.,)°"~" /(E/214.3)"/"—" 


Values of V 7 for the three ratios of 5/t,, are listed be 


low. 
Specimen FPF = 0.366 VY FE 
numbers b/t, (b/t 75 F../F, oie 
1-14 7.3 35.700 0.520 1 aN0 
15-35 14.4 $3,200 0 598 1 908 
40-53 11.3 51,700 0.716 0 9Rn 


Eq. (8) becomes: 

for specimens 1-14 (h, = h — 0.0626) 

m/h,? = 3,570,000/(Bh,/t,?) (Sa 
for specimens 15-35 (h, = h — 0.138b) 

m/h,* = 4,240,000/(Bh,/t,?) Sb 
and for specimens 40-53 (h, = h — 0.1925) 

m/h,? = 5,490,000/(Bh./t,”) Se 
These equations are compared with the test data from 
reference | in Figs. 7, 8, and 9. In the above equations, 
h, is the distance between flange centroids (Appendix 
B). 


DESIGN CURVES 


The development of design curves for predicting the 
ultimate moment of a multiweb box beam in pure bend- 
ing can be simplified by considering a ‘“‘unit’’ of the 
beam as made up of a single interior web and its ac- 
companying effective width of skin (Fig. 2c). If P, is 
the allowable flange load per web, 

P, = F..(2w,t, = 2bdt,,) 
and the allowable moment, .J/, per interior web is 
M = Phy = WwetheFee} [(dtw)/(wels)] + 1f (9) 
Dividing Eq. (9) by hf,” 
MT /(hets?) = 2Fee(We/te) | [(O/te) /(We/ts) te/ts) +14 


and from Eqs. (1) and (5) 


TW RVFGE ea ty) ] 


hd,* jot)" i. Rist)" 


10 


Eq. (10) is plotted for 75S-T6 aluminum alloy (mini 
mum guaranteed properties) in Fig. 10 for those cases 
where f; > ¢,, and 50 > b/t, > 9.7 


7 For the type of construction under consideration, a value ol 
b/t,, of about 9 represents the practical minimum due to clearance 
problems in installing web-flange rivets 
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ULTIMATE STRENGTSA 


A multiweb box having n-webs will develop a bending 
moment equal to n-times the value given by Eq. (10). 
For precise results, however, a correction must be in- 
troduced to account for the fact that the exterior webs 
each have but one effective width of skin. It is there- 
fore necessary to subtract the amount of resisting 
moment due to 2w, of skin. This correction, expressed 
in terms of the parameter .J//h.f,?, is derived in a manner 


ynalogous to Eq. (10), is given by 
Mo /(hes?) = (RV FeyE)/(b/te)" (11) 


ind is plotted for 75S-T6 in Fig. 11. 
As an example of the use of the design charts, assume 
, 79S-T6 multiweb box beam of the following propor- 


tions. 
0.051 in. 
0.125 in 
h = 3.76 m. 
h = 0.74 in. 
nN = & 
} — ) 15 
= 14.5 
h, = 3.66 in. [Eq. (B-1), Appendix B] 


From Fig. 10 17/h,t,? = 1,320,000 Ibs. per sq.in. From 
Fig. 11 1%)/h.t.2 = 1,120,000 Ibs. per sq.in. The allow- 
able bending moment for the box is 


M = (m X 1,320,000 — 1,120,000)/,t," 
= 9,440,000 X 3.66 X (0.125)? 
540,000 in.Ib. 
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APPENDIX A 


The plot of V7 vs. 6/t, (Fig. 6) was obtained as 
follows: For a given material, such as 75S-T6 aluminum 
alloy, values of V7 as a function of stress were ob- 
tained from the stress-strain diagram and Stowell's 
Equation [Eq. (6 Rewriting Eq. (1) with C, = 0.366 
(angle sections supported on two edges), 

Fo/ Fey = 10366/(6/1,)' | V E/Fe (A-1 
The quantity V 7 can also be expressed as a function of 
F../F-, by assuming the stress, F,, obtained from the 
stress-strain diagram, equal to F,, and dividing all 
values of F,, by the yield stress. Thus Eq. (A-1) be- 
comes 


Fic] Fey = (0.366 V E/F.y)/(6/te) “* = f(V) (A-2) 


from which V 7 is expressed in terms of 5/f,, as in Fig. 6. 


APPENDIX B 


The effective height, /,, of the element shown in Fig. 
2c is obtained from the relation 


he =h + 2y (B-1) 


€ 


where 
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_ 6 (3 b hy 1) | bs we (4)")t 
~" 9W2 hn” te /Lte * 4 \ted JS 
(B-2) 


(B-2) was derived assuming the centroid of the 
angle at a distance b/4 from the median line of the leg 
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of the angle that is riveted to the cover sheet. The 
ratio w,/t, is obtained from Eqs. (5) and (1) as 


w,/ts = 1.365k(b, lw) 


where, as before, k 1.03 ws VE/Fey 


Design Charts for Transient Temperature Distribution 


(Continued from page 762) 


the temperature potential for heat-transfer from the 
boundary layer to the wall is also a function of time. 
The design charts represent many possible solutions 
not considered here for lack of space. The generaliza- 
tion of the present analytical solution of the one- 
dimensional heat flow problem to other cases will be 
given in a separate publication, but the present design 
charts will not be changed. 

The methods of using these design charts are illus- 
trated by two detailed numerical examples, based on a 
previous study! of wedge-shaped wing moving at 
supersonic speeds. These examples show that the 
design charts yield answers that are as accurate as those 
obtained previously by laborious numerical integration 


of the heat-flow equations. Furthermore, the time 


required to obtain the answers with the use of the de- 
sign charts is a small fraction of the time needed in 
numerical integrations, since the entire history need 
not be computed when using the design charts. It is 
also evident that the two-dimensional heat flow in a 
thin object moving at supersonic speeds can be solved 
using the design charts by obtaining a series of solu- 
tions for one-dimensional heat flows adjacent to one 
another, each normal to the chord of the object. 


REFERENCES 
1 Kaye, J., 
Flying at Supersonic Speeds, Journal of the Aeronautical Sciences, 
Vol. 17, No. 12, pp. 787-808, December, 1950 
2 Kaye, J., A Table of the First Eleven Repeated Integrals of the 
Error Function, J. Math. and Phys., M.I.T. (to be published 


The Transient Temperature Distribution in a Wing 
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A Rigorous Method for Finding the Lift of a 
Certain Class of Airfoils and Remarks on the 
Meaning of Schrenk’s Approximate Rule 


OTTO LAPORTE* ann HIDEO YOSHIHARAt 


Unwersity of Michigan 


INTRODUCTION 


S IME YEARS AGO, one of the authors! published a 
method for finding on the basis of the lifting line 
theory, the exact lift distribution for a certain rather 
general class of airfoils. With the exception of the 
elliptic plan form this had previously only been possible 
by elaborate numerical methods. While for a given 
plan form these latter methods will still have to be 
used, it seems desirable to have exact solutions avail- 
able for a plan form family characterized by two 
parameters, namely, aspect ratio and taper. 

The plan form which is investigated in the following 


has zero sweep and is given by 
c(y) = @V [1 — (2y/b)?] V [1 — t(2y/d)?] (1) 


In this formula, y is the span coordinate, } is the span, 
cis the chord at y, cy is the maximum chord at the 
center of the wing. The occurrence of a second con- 
stant /, which may be called the taper constant, is new. 
For ¢ = 0 the above plan form degenerates into the 
elliptic one; for ¢ between 0 and 1 wings, which are 
more tapered than an ellipse; and for ¢ between 0 and 
—1 wings blunter than the ellipse result. In the paper 
cited, only untwisted wings were discussed. Since 
that time it has been possible (1) to shorten and sim- 
plify the analysis given there to a certain degree, and 
2) to extend the theory to the case of twisted wings 
i.e., to cases where the angle of attack is a function of 
vy. These developments constitute the first part of 
this paper. 

In this way, exact lift formulas are available for the 
first time for a comparatively wide class of plan forms. 
It seems, therefore, an inviting prospect to investigate 
to what extent the approximate method of Schrenk? 
compares with the exact formulas. Although this 
approximate method is widely used in the industry, 
any theoretical justification has been lacking. 

It is seen that for untwisted wings Schrenk’s rule 
arises out of a simple approximation performed on our 
rigorous formulas. Since the meaning of this approxi- 
mation can be recognized in general, there results a 
generalization of Schrenk’s rule which is independent of 
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plan form and applicable to any kind of twist or even 


discontinuous a distribution. 


DESCRIPTION OF THE METHOD 


According to Trefftz,* the problem of the wing re- 
garded as a lifting line may be reduced to the following 


potential problem in the y, z plane perpendicular to 


the direction of flight: a potential function ¢g(y, 2) 
satisfying Laplace's equation 
Ag = 0 (2) 


which possesses a line of discontinuity or a branch cut 


on the y axis for y between —(b/2) and (6/2). This 
line of discontinuity is the projection of the wing of 
span 6 on to the y, z plane. On the upper side of this 


slit g(x, y) must satisfy the following boundary condi- 


Oy S¢(y, 0) 
02 /,=0 mc(y) 


Furthermore, 


tion :f 


— Va(y) (3) 


—¢(y, —()) (4) 


gly, 


+0) = 


while the normal derivative is continuous. In Eq. 
(3), V is the velocity of the air at infinity, a(y) is the 
angle of attack which may vary arbitrarily as a func- 
tion of span, c(y) is the chord or plan form and m is the 


slope of the lift curve as it results from the infinite span 


theory. The connection between the potential and 
the section lift at y or the downwash velocity is given 
by 

l(y) = dp Vey, Q); wy) = [O¢(y, z) Oz]. 9 (5 


In potential theory a problem such as the present 
one is called a boundary value problem of the ‘‘third 
kind,’’ 


and its normal derivative (in fact, a linear combination 


because a linear combination of the function 


with a variable coefficient) is to assume a given value. 
It may be noted, for future reference, that if the func- 
tion ¢ itself is to assume given values we have a bound- 
ary value problem of the “‘first kind,’’ while if the nor- 
mal derivative is prescribed on the boundary the 
boundary value problem is said to be of the ‘second 
kind.”’ 

The reason why the present problem does not lend 
itself readily to analytic attack is that a variable coef- 


t See e.g., reference 4, pp. 98-112 
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ficient occurs in Eq. (3). Thus in the past, solutions 
have usually been obtained by more or less numerical 
methods. (The only exception to this is the elliptic 
wing.) The idea underlying the present paper is the 
following: if in the y, z plane a conformal transfor- 
mation to new variables, say X and J, is carried out 
the differential equation (2) is not affected. Because 
of the occurrence of the first derivative, however, the 
boundary condition (3) will change; in fact each trans- 
formation (y, 3) > CY, V) will cause a new variable co- 
efficient to appear both in the ¢ term and on the right- 
hand side of Eq. (3). We propose to carry out a con- 
formal transformation and then to choose such a chord 
function c(y) that in the new variables a boundary con- 
dition with constant coefficients results. This new 
problem will then be treated with standard methods. 


First TYPE OF CONFORMAL MAPPING OF THE TREFFTZ 
PLANE INTO A RECTANGLE 


We regard y and s as one complex variable by put- 


ting? 

u=y+2 
It should be remembered that the trace of the wing 
is a slit on the real axis extending from y = —(b/2) to 
+(b/°2). The region outside this slit will now be 


mapped onto a complex plane 
Z=X+1Y 
by means of the function* 
u = (b/2K)ns[(2/r)KZ, K] (6) 


In this equation Z is the complex variable. & is a con- 
stant parameter called the modulus which we will as- 
sume to be real and between 0 and 1, A is the complete 
elliptic integral of the first kind of k. ns(x, A) is the 
reciprocal of the elliptic function® s(x, k) of argument 
x and modulusk. The mapping of Eq. (6) is illustrated 
in Fig. 1 by denoting corresponding points in the two 
planes by identical numbers. The essential features 
arise as follows: 

(a) Z = O gives, since sv of vanishing argument 
vanishes, the infinitely distant point w = ©, 





* From here on the presentation is different from and simpler 
than the original derivations published in Quart. Appl. Math., 
1944. 
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(b) The upper or lower sides of the rectangle Z = 
X + 1(7/2) (K’'/A), where AK’ = A(k’) with &? + 
k’? = 1, map, respectively, into the lower or upper 
side of the slit along the real axis of the u plane between 


y= +(b/2). Thisis seen by calculating 


u = (b/2k)ns(2/m)K[X + i(e 2) (K’'K 


le(m 2K)] 


A small increment ¢ is added here, to enable us to in- 
vestigate not only points on, but also points near, the 
sides of the rectangle. We obtain® 


u = (b/2) [sn(2/m)KX F 


tecn(2/rm) AX dn(2 rm)KX 7 


Thus both edges of the rectangle become the real ) 
axis from —(b/2) to +(6/2). We have therefore to 
imagine that the w plane is cut open along this piece of 
its real axis. 

(c) We turn to the sides of the rectangle by calcu- 
lating the « values which go with Z = +(m 2) + 7¥: 


u = (b/2k)ms[tA + 1(2/m)AY] 


that is to say, 


y = +(b/2k)dn[(2/m)AKY, K’] 7’ 
Since dn of a real argument is real (and even) u will be 
real. The corners of the rectangle ) = +(7? 
(K’/K) must give the same map points as under (b); 
since dn(+K’, k’) = k we get, indeed, 1 = +)? 


As Y diminishes from + 7A ’/2K to zero, dn increases 
from k to 1. The points in the w plane which corre- 
spond to the midpoints of the two vertical sides of the 
rectangle are therefore the points 1 = y = +0 2k. 
Since & is less than 1, they lie beyond +62. One 
sees in this way that the branch cut on the real axis has 


to be lengthened to +(b 2k). 


SOLUTION OF THE PROBLEM IN THE Z PLANE FOR 
TAPERED WINGS 


The boundary value problem in the ‘‘physical’’ w plane 
was stated 1n the first section. Since by Eqs. (6) this 
plane is mapped into the Z plane, the boundary value 
problem should be restated in terms of this new variable. 
It is clear that the solution of the problem requires 
the construction of a velocity potential g(X¥, J 
satisfying Laplace’s equation (with Y and J’ as inde- 
pendent variables) which fulfills different boundary 
conditions on the various sides of the rectangle. 

Eq. (3) states the boundary condition for z = +0, 
y| < (6 2). Since this line becomes the lower side 
of the rectangle, Y = —(7/2) (K’/K); |X| < (1/2), 
it is clear that now the boundary condition is 


(5 4 S ( ) (x wh ) 
oY } rK’/2K MC dZ } eK’ oe sided 2K 
r (“*) 
—Va : >) 
a2) y = — 2K'/2K 


with 
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FINDING THE Lift OF A 


du bd : ae 
: = - 92 — KX 9 
dZ } eK’'/2K 2 dX us 
[he latter results from Eq. (7). In order that this 


boundary condition have constant coefficients we 
choose c(y) proportional to du/dZ. To see what plan 
form this leads to it becomes necessary to write Eq. 
9) in terms of the original span coordinate y. This 


is readily done using Eq. (7) and gives 
ju \ hK Dy\2 
7) V1 - =} x 
dZ/ +K'/2K 7 5 


\1 - HO) 9’) 


[hus we are led to Eq. (1), but since k is supposed to 
be real and a genuine fraction, only wings of taper 
sreater than elliptic result. Calling c) the maximum 
chord at the center, the area S of the wing is calculated 
readily* in terms of the two complete elliptic integrals 


Eand A: 
5 = boG(k) 
; ereee 10) 
Gk) = (1 3k*) [1 + RE — (1 — RAN 
The aspect ratio @ = b°/S and the average chord ¢ 
become, with this notation, 
Q@ = b/(Goo): c= Goo (10’ 


Substituting Eqs. (9) and (9’) into Eq. (S) and using 
Eq. (10’), the boundary condition at the lower side of 


the rectangle of Fig. 1 becomes 


O¢g SQ bd z , 
ci~- tens Va -—w- KX (1 
0} 7m 2 dX 1 


\t the upper side the sign of the second term has to be 
reversed in view of Eq. (4). From here the solution of 
the problem proceeds according to standard methods. 
Let the right-hand side of Eq. (11) be expanded in a 


Fourier series for Y| Z (2/2): 


«me KX = 5 a, cos (Oe + DX 12) 
dX T n=0 
[he assumption will be made here that the angle of 
attack is a symmetric function of y and therefore of 
\, otherwise sine terms appear in the series. Unsym- 
metrical a distributions will be reported on briefly 
later. 
Phe following assumption as to the form of ¢ is now 
made 
(AX, V) = So d,sh(2n + 1)¥ cos (Qn + 1)X (13) 
0 
Because of the fact that ¢ has to satisfy Laplace’s 
equation the dependence upon VY has to be that of a 
hyperbolic function. The hyperbolic cosine has to 
be ruled out because of Eq. (4). After substituting 
Eqs. (13) and (12) into the condition (11) for Y = 


* Using, for example, Hancock, p. 62, second formula from the 


bottom 


-BRTAIN CLASS OF 


AIRFOTILS 7SY 


—nrk’ 2A the following value for the coefficient 


results: 
l I "ha 
2 (2n l D 
rh’ 
D, = ch(2n + 1 : 14 
ae 
S@Q GK rh’ 
sh (2n + | 
rm2n + 1 2A 


Eqs. (13) and (14) contain the complete solution of the 
problem. The lift on the airfoil itself is obtained from 


Eq. (13) by putting } tA’ 2K. If a sectional lift 
coefficient c, is defined by 
tpl¢ 1 (2) pV°c LS 
then c,is found to be 
~ a rh’ 
ri 1a > h(2n + = 
0(27 + 1)D 2A 
cos (2m + 1) i 


To this formula we add the relation between _Y and the 


span coordinate y: 


y = (b/2)sn(2/mr)KX 1G 
which arises from Eq. (7 
The downwashw = (0g, 0s),—» can be obtained from 


Eq. (3) by substituting Eq. (15) into it. 


VARIOUS CASES OF TWIST 


The coefficient a, in Eq. (15) remains to be deter- 
mined. A twist of the form 


will be assumed. For untwisted wings only the first 
term is used; otherwise 6 and y may be adjusted to 
fit a given structural or aerodynamic twist. This ex- 
pression has to be written as a function of X which 
can be done with Eq. (16). The left-hand side of 


Eq. (12) may thus be put into the following form: 


d 2 d 2 
a _ SH KX — ~<«s (« 7 & SH AN 
dX 7 dx T 


To obtain the coefficients a, of Eq. (12) one needs there- 
fore the Fourier expansions of the odd powers of the 
sn function. These are standard developments.’ The 


results are 


, ae T sin (2m + 1)X 
ss - KX = — Z, —— 
T Kkn=0 sh(2n+1)(7 2)(K’ A 
») av ») r 
rT |. sin (27 + 1). 
s*- KX = — 2, £,- IS 
T Kkn=0 sh(2n+1)(r 2)(K' A 
ee rm. sin (2m + 1)X 
seé- KX = —2Z, F, = 
T Kkn=0 sh (2n+1)(9 2)(K’. A) 
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2 
} Y 
(3) 
f ‘ (6, ln (8) 
Bi | 
k T 
: 
ae )} | (2 + ‘ | ” 
” . ~ t y — — —+-. — — >X 
(4) } 2) 
, }@ 
¢ b - * 
en Mi IR—-| 
u=y +l2 PLANE Z=X+i¥ PLANE 


Fic. 2 


E, = (1/2k?){1 + k? — [(2n + 1) (9/2K)]?} 

= (1/4k*)}(1/2)(3 + 2k? + 3k4) — 
(5/3) (1 + R?) [(2n + 1) (/2K)]?2 + 
(1/6) [(2n + 1) (4/2K)]*} 


| 
| 


(19) 


The coefficient a, is obtained from this, and putting it 
into Eq. (15) the lift coefficient becomes 


dr — | 
= i¢ 0 _ re ‘OS ? E 3 y 
C Kk t dep (a + BE, + yF,,) cos (2n+1)X (20) 


The convergence of this series is quite rapid because 
the denominator D,, of Eq. (14) increases exponentially 
with m. For untwisted wings 8 and y should be put 


equal to zero. 


SECOND TYPE OF CONFORMAL MAPPING 


Consider now, in place of Eq. (6) the function 


5) 
u = (b/2k) ds (7x2, #') (21) 
us 
where k’ = (1 —&? ) is the complementary modulus 


which is also real and between 0 and 1. K is, as be- 
fore, the elliptic integral belonging to k, not to k’. 
The function (21) maps the u and Z planes upon one 
another in a fashion illustrated by Fig. 2. Again 
corresponding points in the two planes are designated 
by identical numbers. In the present case the u = 
y + zz plane must be cut by a cross shaped branch cut 
which possesses the points y = +6/2 and z = +(b/2) 
(k’ k) as end points. This is readily seen by two sub- 
stitutions into Eq. (21). 

(a) The two vertical sides of the rectangle corre- 
spond to the two horizontal arms of the cross. For 
when putting Z = +(7/2) (K’/K) + 7¥ we get* 


u = +(b/2)cn[(2/mr)KY, k] (22) 


u is therefore real. In particular for Y = 0: y = 6/2 
and Y = +2/2: y= 0. 

(b) The two horizontal sides of the rectangle map 
into the vertical sides of the cross. Put Z = +7(2/2) 
+ X into Eq. (21). The result ist 

* Hancock, p. 33, second set. Primed and unprimed quanti- 


ties must be interchanged. 
+ Ibid., p. 33, last set of formulas on page. 
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OR se 
=> 9 n ( KX. r’) 


cd 
2k T 


u is imaginary, and the sign shows that the upper side 
of the rectangle gives the bottom branch of the cross 
and vice versa. The length of the vertical arms, (4/2 
(k’/k) may be greater or less than that of the horizontal 
arms accordingly ask’ 2 k. It will be seen that only 
one of these cases, that of k’ > k, has practical intercst, 


SOLUTION OF THE BOUNDARY VALUE PROBLEM FOR 
BLUNT WINGS 


Boundary condition (3) on the lower side of the wing 
slit becomes, when transformed into the Z plane with 
Eq. (21), 


(oe) 8 (“*) (= ')- 
OX/ x= 5K'/2K ’ me AZ] x = «K'/2K *\ OK’ ~ 


or (“) 
1Va : (23 
dZ/ x = 5K'/2K 


As previously, the chord function c(y) is chosen pro- 


portional to du/dZ which now is, using Eq. (22), 


du bd we 
2 =-—-1 -cn KY,k 
dZ X=nrK’'/2K 2 d) us 


_ ,bOK bali (2) ; 4 k? (7 : 
is T \ b \ eo 


We therefore choose as chord 


Av) = ¢ l — (7) l + k* (7”) ») 
i \ b \ k’*\ b - 


which gives a wing blunter than the elliptic wing. 
Since for k > k’ this function has a maximum at a 


place between the center y = 0 and the tips y = +) 2 
we shall only consider moduli k for which 
G< £ <2 
The entire family of plan forms—i.e., both the 
tapered ones given by Eq. (9’) and the blunt ones 
given by Eq. (24)—are shown for the arbitrarily selected 
aspect ratio of 8 in Fig. 3. 
We also note that for .S, @, and @ we have now 


S = beooH(k)/k’'; @ = bk'/(H(k)col; 
€ = H(k)ao/k’ (25 
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Fic. 3. Plan forms of @ = 8 for various K values. 
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FINDING THE LIFT oF A 


where //(k) is another standard elliptic integral 
H(k) = (1/3Rk?) [1 — R2)K — (1 — 2k?)E] 


For the limiting value k = 2 ? which interests us 
most A(k) has the numerical value 0.61802. The 
wing shape for this case is the full line in Fig. 3. 

In accordance with the fact that now the condition 
23) has to be satisfied on the two vertical boundaries 
of the rectangle, the following hypothesis as to the 
form of ¢ is made: 

¢g = b,’ch(2n + 1)X sin (2n + 1)¥V (26) 
n 0 


From then on the calculation runs quite analogously to 
that of Eq. (5). The final formula for the section lift 
coefficient, in case of a twist given by Eq. (17) is 
= ai » (ao + BE,’ + yF,') X 
Kk n=0D,,' es 
sin (2n + 1)¥ (27) 


The relation between Y and the span coordinate y is 


b 2 
é TT 


The following abbreviations, analogous to Eqs. (14) 
and (19), are employed :* 


D,’ = sh(2n + 1) SN 
2K 
8 @ HK ch(Qn + 1 7K! 
am2n + | 2K 
EE’ = J_j + 2k? + | (on + 1) =|" (29) 
2k? | 2K IS 


L jl 5 
Bt = <— (3 — Sp? + SR) — 1 — 2k?) > 


2 | 4 
(2n + 1) + (2n + 1) af l 
2K 6 2K | 4 


The convergence is as rapid as that of Eq. (20). For 


« 


untwisted wings 8 and y vanish. For twisted wings 
the coefficients -,’ and F,,’ simplify somewhat if one 
restricts oneself to the interesting limiting case k = 
») 


REMARKS ON DISCONTINUITIES OF THE VELOCITY 
POTENTIAL ¢ IN THE PLANE OF THE ORIGINAL 
VARIABLES y AND 2 

As stated in Eq. (4) it is desired that ¢ have a dis- 
continuity along the y axis when y is between +0/2. 
(The normal derivative, however, should be con- 
It is immediately seen that our ¢ solutions 
But the con- 


tinuous. ) 
(13) and (26) behave in this fashion. 
formal representations (6) and (21) introduce, as Figs. 
1 and 2 show, further lines of ramification, namely, 

* The coefficients E,’ and F,,’ are essentially the expansion 
coefficients in the Fourier series for the cube and the fifth power 
of cnu 


CERTAIN CLASS OF 
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to +6, 2k or, in the, 
It should be 
verified that on these ‘‘unphysical’’ cuts the potential 


either a branch cut from +) 2 
blunt case, from +0 to +(b/2) (k’/k). 


yg or its derivatives do not suffer discontinuities. 

For tapered wings, Fig. 1 indicates that opposite 
points on, say, the cut between 6/2 and 6/2k corre- 
spond to points whose Z is 7/2 + ty. Eq. (13) shows 
that for such values ¢ vanishes. Turning to the de- 
rivative we see that, by virtue of the Cauchy-Riemann 
equations 


(sx )5 ¥) == (2) © (o5)(3-») 


Since the left-hand side is odd in Y, and the second 
member of the right is, by virtue of Eq. (7’), likewise 
odd in Y, the first member must therefore be even in 
Y—i.e., continuous for points on opposite sides of the 
cut 6/2 < y < 6/2k. 
Laplace’s equation and the Cauchy-Riemann equations 


Similarly one proves, using both 


that all even derivatives of g with respect to z vanish 
on the cut, while all odd ones are equal for opposite 
points. The proof for the vertical cuts of the “blunt”’ 
mapping is analogous. 

Thus it is seen that, although the mapping possesses 
branch cuts in addition to the physical one —(b/2) < 
y < (6/2), the velocity potential ¢(y, z) has the re- 
quired discontinuity only at the latter cross cut while 
behaving perfectly continuously at the former ones. 


FURTHER RIGOROUS CALCULATIONS 


For the sake of completeness, two extensions of the 
calculations presented in the previous sections should 
be pointed out. 

(1) The first of these is concerned with the limiting 
case of extreme taper when the plan form is given by 


c(y) = afl — (2y/d)?] 


This case is treated in detail in Section 8 of reference 
1. It is of some interest, for it contains the rigorous 
solution of the lift problem for a wing whose tips are 
not rounded off, but rather pointed. The velocity 
potential g and the lift coefficient now appear in the 
form of series proceeding according to the roots o,, of 
the transcendental equation t 


o, COS o, + asino, = 0; a = (47/3) (@/m) 


The result for the lift coefficient c, of an untwisted 
wing turns out to be 


° 9 ¢g, 
= at+1 Onl@*® + o,*) “*/1 — 2y/b 
C: = 8a@ > (—1) 
n=1 


a(a+1)+ | 1+ 2y b 


The roots o, can be obtained graphically with ease. 
(2) In the doctoral dissertation of the junior author® 
still another case of twist is calculated, namely, one 
where the angle of attack is a linear function of span. 
(3) Compact expressions for overall quantities, such 
as total lift, total drag, and rolling moment, have been 


+ The above a value should be substituted for Eq. (51”), p 
244 of reference 1. 
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obtained, and their dependence upon taper, aspect 
ratio and degree of twist calculated. For this too, 
see reference 8. 


(4) Tables and graphs of the dependence of the 


sectional lift coefficients upon span may also be found 
there. 


AN APPROXIMATION LEADING TO SCHRENK’S FORMULA 


As stated in the introduction, the rigorous formulas 
obtained will be approximated in such a fashion that 
Schrenk’s rule results. For the sake of briefness this 
process will be explained using only the formulas of 
the tapered wing. The formulas for the blunt wing 
can be dealt with in an analogous fashion. We shall 
also restrict ourselves to untwisted wings, putting 8 and 
y in Eq. (20) equal to zero. 

The convergence of the series (20) depends essentially 
upon the denominator D, which, from its definition 
(14), is seen to consist of two terms, both rapidly in- 
creasing functions of 27 + 1. We propose to replace 
the denominator D, of Eq. (20) by two simpler ex- 
pressions. We shall 

(a) Replace ch(2n + 1) (m/2) (K’/K) by A(2Qn + 
1)—!sh(2n + 1) (2/2) (K’/K). 

(b) Replace (2m + 1)sh(2n + 1) (m/2) (K’/K) by 
\—! ch (2n + 1) (7/2)(K'/KR). 

For small values of this is a perfectly plausible 
approximation; however, for large values of n, while 
both expressions increase exponentially they deviate 
from one another more and more. As far as large n 
values go, it would seem more plausible simply to re- 
place the ch by the sh function, and vice versa. If, 
nevertheless, the above method of approximation is 
chosen, it is in view of the fact that terms with small 
n contribute overwhelmingly more in Eq. (20). The 
numerical constant A will be determined later so as to 
give the best approximation in the mean. 

When replacing, according to (a), the rigorous ex- 
pression for the lift per unit span / is changed into an 
approximation which will be denoted by /;. This be- 
comes, according to Eqs. (15) and (20), 


2a app Vb 
i= —— . — 2) (2n +1) X 
KkX + (8/r) (@/m)GK n=0 
cos (2n + 1)X 
’ tae (30 
sh(2n + 1) (w/2) (A’/K 
This is, according to the first formula (18), proportional 
to 


cn(2/m)KX dn(2/r)KX 
and therefore, according to Eqs. (9), (9") and (16), 


1A ap Vb 
= ccc: Oe 
xt X\ + (8/7) (@/m)GA 


\ [1 — (2y opal — HG) | (30’ 


On the other hand, replacing according to (b) and call- 
ing the resulting approximate lift /2, 

27 app V*br 

KkX + (8/2) (Q/m)GK 
= cos (2n + 1)X 


» 
n=0 ch (2n + 1) (7/2) (KK 


9) 
(31 


The series is the well known expansion® of cn(2/7)AX. 
Therefore, 


2rayp V7b Dy\: 
I = RVI - (= (31’ 
A + (8/7) (@/m)GK bh 


To fix XX most advantageously, it seems reasonable to 
calculate the total lift by integrating Eqs. (30) and 
(51) along the span, and to equate the obtained ex- 
pressions. We get first from Eq. (30) 
4K ay pV 


L a = - § 
aX + (8/7) (@/m)GA 


where G is the elliptic integral introduced in Eq. (10). 
Then from Eq. (31) 


2raop Vb ,* 


2 = = 
A + (8/r) (@Q/m)GA 4 
Equating these two expressions we obtain for X, 
\ = 8KG/r° (32 


while /; and /, assume the following form: 


l . Tayh c(y) 
rf = pl . “2 
2 1+ 7rQ/m c¢ 
(33 
l FE Tad V1 — (2v/b)? 
ls = pl >i = 
a 1 + 7rQ/m w4 


If we now take as a less biased approximation the 


arithmetic mean of the two, 


we obtain exactly Schrenk’s formula for untwisted 
wings.* For blunt wings the same formula—but with 


* See first formula on p. 3 of Schrenk’s paper quoted in the 
introduction. A correction constant A, used by him to adjust 
the dependence upon aspect ratio turns out to be zero in our cas¢ 
when J is fixed as in Eq. (32). 
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FINDING THE LIFT OF A 


the chord of Eq. (24)—results, after replacements 
corresponding to (a) and (b) have been made in the 
denominator D,,’ of Eq. (29). 

In Figs. 4 and 5 the coefficients 1/D, of the exact 
solution together with the Fourier coefficients of the 
series (30), (31) for /, and /, are shown for two values of 
the modulus k. 


GENERALIZATION OF THE PROCESS EMPLOYED ABOVE 
ro ARBITRARY PLAN FORMS AND TWISTS 


When looking at the rigorous solution (19) we see that 
its coefficients E,, F,, and D, owe their appearance 
to different causes: /, is due to the quadratic twist, 
F,, to the quartic twist, while D,, arises from the funda- 
mental boundary condition (3). It is this latter co- 
efficient which was replaced either by (a constant 
times) its second term or by (a constant times) its 
first term. But since the first term of D, is due to the 
first term of the boundary condition (3), and the 
second term of D, to the second term of Eq. (5), it 
may be said also that in the preceding paragraph we 
have replaced the rigorous boundary value problem 
(2), (3), and (4) by the following problems: 

a) The boundary value problem of the first kind 
with the boundary condition for —b/2 < y < b/2 

Vam ; 
gi(yv, 0) = - c(y) 
ii Al yay 

b) The boundary value problem of the second kind 

with the condition for —6/2 < y < b/2 


(=) nN Va 
Os / s=0 1+) 
In addition, both g; and g» have to satisfy Eq. (4) and 
have to be solutions of Laplace’s equation (2). The 


constant is fixed by the requirement that the integrals 
along the span of ¢; and of g: agree. The Schrenk 
approximation is the arithmetic mean of ¢, and ¢». 

The above problems may readily be solved for any 
chord or angle of attack distribution. But inasmuch 
as we are interested only in lift—i.e., in ¢ (y, 0)—the 
value of ¢;(y, 0) is already fixed by the above condition 
(a). Problem (b) is seen to be independent of the 
plan form and can therefore be solved as if an elliptic 
wing were being studied. 

We shall illustrate this for the case of an untwisted 
wing by rederiving—for an arbitrary plan form—the 
two approximate lifts (33) whose arithmetic mean 
Schrenk’s formula (34) employs. /; is taken immedi- 


itely from (a) and gives 


1 2 am , - 
1, = ~ pV? —— c(y) 39) 
2 ae 
and a total lift 
I l me am 
Yes 5 p re . 


¢: and /, we anticipate to be elliptic since the down- 
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Lot K= YG 
— 
aa, st" m7 20+ orf 
om /6+ tTAle A GR COSEL 
+ 
~~ Dn 
N 
+ 24) "sech 
4 ™ ee 
.~ 
O04 4 > 
= Ee O 2 4 6 8 < 
/7 
Fic. 5 
wash derivative is constant. The conformal trans- 
formation* u = (b/2)ché maps the exterior of the wing 


slit into the semi-infinite strip bounded by the points 
o, 0, 271, ~ + 277. The two sides of the cut from 
—b/2 to b/2 become the imaginary axis from 0 to 27i. 
Boundary condition (b) becomes, in terms of & and 


n l(t = —&+ In), 


(=) l"aBr , 
= - Sin 
dt /.—0 21 +2 


$ 


The solution of this boundary value problem is 


l"abv : a 
go = sin 7e 
2(1 + A) 


‘ 


which gives as lift 


2p V%abd Dy\2 
) = - : 30 
1+A \ (= ) 


and as total lift 
l 172 mab") 
2 1+A 


Equating the Z we find 
\ = (m/z) (S/67) = m/rQ@ 


After substituting this value for \ into Eqs. (35), (36), 
these / values are seen to agree with Eq.’ (33). 

A last remark on discontinuously varying a distri- 
butions, which occur in the case of flaps or ailerons, may 
be added. The above cited paperst have several fig- 
ures which demonstrate that Schrenk’s original recipe 
gives rise to discontinuous lift distributions. These 
being unreasonable, more or less arbitrary methods are 
given to ‘‘fair out’’ the discontinuities. Let us see 
what our approximations (a) and (b) give instead. The 
function g; being proportional to a is also discontinu- 
ous. But ge will be perfectly continuous even though 
its normal derivative may jump. The arithmetic 
mean of g; and ¢g» will thus show only a discontinuity 


of half the amount. 


* See reference 1, Section 2, p. 234 


t See reference 2 
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Remarks on the Velocity of Sound 


Max M Munk 
The Catholic University of America, Department of Aeronautical 
Engineering, Washington, D.C. 


May 15, 1955 


HERE IS A NEED for a wider spread understanding of the 
y jarerth meaning and definition of the velocity of sound 
This is a fundamental quantity in the theory of gas dynamics 
Whoever wants to understand the theory has first to understand 
the velocity of sound 

We learn to understand and to compute the velocity of sound 
by scrutinizing how it is first introduced into the equations 
The initial and most basic equations of gas motion do not con- 
tain the velocity of sound. It is introduced as a later step to 
get the theory under better control. 

The theory deals with the motion of a continuous fluid in 
three dimensions, not necessarily a steady motion. Shear forces 
are assumed to be absent. This is the only limiting assumption. 
Viscosity is excluded only as far as it would introduce shear 
Heat conduction and addition or removal of heat are 
That by itself does not lead to shear. The 


forces 
not excluded at all. 
gas is not necessarily perfect, not barotropic, not even homo- 
geneous 

The broad class of shearless fluid motions is governed by the 


system of equations 


+ (grad p)/p = O| 
4 { } 
o + divv = of 


The first of these two equations is a vector equation and is equiva 
lent to three ordinary equations. It describes the dynamic rela 
tion between the velocity v, the pressure p, and the density p 
The second equation expresses the condition of continuity 

Eqs. (1) contain five dependent variables. The system (1) is 
incomplete. It is necessary to introduce a combination of the 
dependent variables as a new independent variable. It is ex- 
pedient and customary to choose the velocity of sound for that 
purpose 

In Eqs. (1) the dot represents time differentiation with La 
grangian independent variables. That is to say it designates 
the rate of change of the quantity in question with reference to 
itself 


Such rate of change has a much more vivid physical significance 


each separate fluid particle preserving its identity by 


than items merely describing the variation of the pressure field 
or density field. 

The square of the velocity of sound, c’, is introduced by way 
of eliminating the derivative of the density, by the equation 


p = p/c? (2) 


The quantity c? relates accordingly to one specific gas particle 
and its doing. It does not relate at all to the adjacent particles 
or to the field distributions and their variations. 
1/c? is a kind of de facto specific volumetric elasticity modulus 
of the particle, a kind of effective hardness of the gas. It is the 
intensity by which it appears to resist the change of its volume 


The expression 


under the action of the pressure, as if all volume changes were 


attributable to the pressure and none to the change of temper 
ature. That is to say, the pressure changes are made responsible 
for the temperature changes and its effects too, whether they art 
or not 


It is often stated in the literature 

c? = dp/dp 3 
Eq. (3) is consistent with Eq. (2) provided the density through 
out the flow is a function of the pressure only. If it were not 
the ordinary differential quotient in Eq. (3) would not even have 
any meaning. If the pressure and the density vary independ- 
ently of each other, only the partial derivative 0p/0p would 
come in the question. That partial derivative has not a unique 
value for any point or particle. In a steady flow the partial 
derivative 0p/Op depends on the direction of the line along which 
in succession the values of the pressure and of the density are 
studied, and their changes compared with each other No such 
lack of uniqueness attends the velocity of sound. Only one of 
the many values of 0p/Op can therefore be equal to the square 
of the velocity of sound. It is computed in a steady flow in the 
direction of the stream line, for the same particle occupies then, 

in succession, all points of that line and of that line only 
Acoustical waves are assumed to be associated with isentropic 
expansions and contractions. Therefore the square of the acous- 
tical velocity of sound is equal to the thermodynamic partial 
derivative 0p/Op for constant entropy. The velocity of sound 
appearing in the equations of the gas motion, including the forms 
of these equations called ‘‘characteristic’’ is not necessarily equal 
is then some- 


to the isentropic derivative. ‘‘Velocity of sound” 


what of a misnomer. Superposing acoustic waves to a gas 
motion changes the flow and its gas properties, and may change 
a steady flow into an unsteady one. That would change the 
magnitude of the applicable velocity of sound, too 

An incompressible fluid having a density variable throughout 
space exhibits many different magnitudes of 0p/Op but its ve- 
locity of sound is, of course, infinite. That illustrates the lack 


of identity between c? and Op /Op. 


Strength Under Combined Tension and 
Bending in the Plastic Range 


J. M. Frankland and R. E. Roach 

Staff Assistant, Structures Technical Group, and Lead Structures 
Design Engineer, Respectively, Chance Vought Aircraft, 
Incorporated, Dallas, Tex. 

May 18, 1955 


NTHONY J. BARRETT has recently published! in the Readers’ 

Forum some interesting comments on our paper? on plastic 
effects under combined tension and kending. He points out, 
that the effects of deformation of cross section 


This circum- 


quite correctly, 
are small and could well have been neglected 
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Fic. 1. Loading history of a curved bar under compression 
loading. 























stance, however, was not self-evident and called for confirmation. 
A statement regarding the negligible effects of change in section 
shape was inadvertently omitted from the published paper. 

We are particularly glad that Barrett has brought up this 
point, since otherwise it might appear that the published curves 
apply only to a material with the stated value of yield strain. 
The complications due to section deformation can be eliminated 
from the expressions in the original paper by setting ey = O in 
Eqs. (3) through (8), (12), (16), and (17). 

Regarding the use of apparent stress rather than true stress, 
there seems to us no advantage to the use of apparent stress. 
This choice enters only into the determination of n, the constant 
characterizing the shape of the stress-strain curve. There is 
little work required to convert the conventional stress-strain 
curve to the true stress-strain curve for the values of strain which 
are of interest, since these are less than the strain required to 
reach maximum load in the tensile test. 

Doubts were expressed in the original paper regarding the ap- 
plicability of our results to the calculation of curved bars under 
axial compression in view of the stress reversals which would 
take place. Barrett, however, points out that this is too pessi 
mistic a conclusion for a material which can be considered elastic 
up to three-quarters of the yield strength, which is true of most 
of the high-strength alloys in use today. He proposes that 
our methods apply for eccentricities larger than one-fortieth of 
the beam depth. He justifies this by consideration of extreme 
fiber stresses. 

It is our belief, however, that the strain history of each point in 
the cross section must be considered, since stress reversals may 
take place in the interior of the beam as well as at the surface. 
Only in this way can one be certain that variations in conditions, 


such as different shapes of the stress-strain curve, and differen; 
redundant combinations of axial load and bending, such as 
eccentric loading of an initially straight member, have not made 
the calculation invalid. 

Considering the case of a curved bar under axial compression 
with an initial eccentricity equal to one-fortieth of the beam 
depth ( Barrett’s approximate limiting case in compression ) and as 
suming the same proportions used in the example in our origina! 


paper,” except for the eccentricity, 


and ey = 6.3 X 10 


For compression loading, Eq. (23)? becomes 


M | (= 3eyl? ) P 
My \h xh?) P, 


For the specific example this reduces to 
M ‘My = (0.15 + 0.5 «) (P Py 


The resulting load history through maximum load for a material 
with a stress-strain curve characterized by m = 21 is shown in 
Fig. | 

From the intersection of this loading curve with the lines of 
constant elongation, m, the strain in the most highly stressed 
fiber, ¢;, can be determined at a number of load levels. The 
values of « can be obtained either by interpolation between 
« curves or by calculation from the slope of the line connecting 
the origin and the intersection of the loading curve and the 
appropriate m-line. That is, 


(4 My 5) 
k=2 — 0.15 
P/Py 


From the value of «x, the strain in the least strained fiber is, 
from Eq. (20)? 


The strains over the depth for the example case are shown in 
Fig. 2 for several values of extreme fiber strain 

Accepting Barrett’s suggestion that strain reversals below three 
quarters of the yield strength take place elastically, and there- 
fore reversibly, it will be seen that up through maximum load 
all strain reversals occur in the elastic core of the section. In 
this particular case, then, the method of our paper ts applicable 
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Fic. 2. Strain distributions over the depth of a curved bar under 
compression loading 





wave 
rotat 
tube 
consi 
by 
refer 
TI 
expe 
inter 


plate 


flecte 
mitte 


phen 








ifferent 
uch as 


it made 


ITeSSion 
beam 
ind is 


original 


iterial 
wn in 


nes of 
ressed 

The 
tween 
-cting 
d the 


er is, 


vn in 


three 
here- 
load 
In 


‘able 


der 





TABLE | 





Values of Stress and Total Strain at Which the Plastic Strain is one 
per cent of the Elastic Component of Yield Strain 
n 7 1] 15 21 
J; 0.5 0.711 0.778 0.836 





in compression certainly down to an eccentricity of one-fortieth 
the beam depth, in confirmation of Barrett’s remarks 
[he level of stress below which stress reversals are permitted 


iould more reasonably be stated in terms of an acceptable plastic 


strain. With this criterion, the value of allowable stress (or 
strain) for stress reversals will depend on n If it is assumed 
that a reversal of plastic strain no greater than one per cent of 
the elastic component of the yield strain should occur, then a 
value of the plastic component of strain, m, equal to 0.01 would 
define the limiting stress. The limiting values of stress for this 
assumption are shown in Table 1 for several shapes of the stress 


strain curve 
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Some Comments on “‘A Theoretical and 
Experimental Study of Shock Tube Flows’’* 


Darshan S. Dosanjh 
Research Associate, Institute for Fluid Dynamics and Applied 
Mathematics, University of Maryland, College Park, Md 


May 21,1955 


ECENTLY GLASS AND PATTERSON! have published some posi- 
R tion vs. time schlieren records of the interaction of a shock 
wave with screens. These were recorded by the technique of a 
rotating drum camera in conjunction with a wave interaction 
tube. They observed that after the interaction, the wave system 
consists of a reflected and transmitted shock wave separated 
by a contact surface and a vortex (see plates 18-19 in 
reference | 

The author? 3 also conducted, independently, some similar 
experiments in a shock tube to study the various aspects of the 
interaction of traveling shock waves with grids (perforated 
plates, Fig. 1). A combination of optical (shadowgraphic) and 
hot wire techniques was used for experimental recording. The 
wave model for shock wave — grid interaction is reproduced in 
Fig. 2. The grids were mounted between two glass sections 
about 16 ft. downstream from the diaphragm location. By this 
arrangement, it was possible to simultaneously record the re 
flected shock wave upstream of the grid, as well as the trans 
mitted shock wave, turbulent contact surface and any allied 
phenomena downstream of the grid. (For details see references 
2and 3 

The presence of a vortex as reported by Glass and Patterson! 
is not substantiated by our experimental observations. A large 
number of shadowgraphs of the shock wave — grid interaction 
were taken systematically for various incident shock strengths 
at controlled delay time settings. No vortex, of the nature re- 
ported in reference 1, was observed. What was, however, con- 


Ihe experimental work on which this article is based?) ? was carried out 

the author during 1953 at the Department of Aeronautics, Johns Hopkins 

University, Baltimore, Maryland The work was sponsored by the United 
States Air Force, under Contract AF18(600)-757. 

The author gratefully acknowledges the helpful discussions with Francis 

H. Clauser. Richard Swartley’s help with the experimental work is appre 








Fic. | The grid used for the investigation of shock wave-grid 


interaction. 1/4-in. thick perforated steel plate with reamed 


5 in. diameter. Solidity ratios = 0.25 





holes of 0.2 


sistently observed was that for incident shock waves of strength 
S > 3.0 where S = the ratio of the pressure behind the incident 
shock wave to the pressure ahead of the incident shock wave, 
there emerged, following the transmitted shock wave and turbu 
lent contact surface, an auxiliary shock pattern out of the grid 
(Fig. 3a Notice that the reflected shock wave and this auxiliary 
shock pattern are similarly faced and the transmitted shock wave 
is oppositely faced to them. 

A controlled set of shadowgraphs, with longer time delays, 
established that this auxiliary shock front was being washed 
downstream at very slow speeds (about 1/2 in. in 1.3 millisec., 
Fig. 3b. The upstream face of the grid corresponds to zero de- 
lay in both of the shadowgraphs reproduced here It is felt 
that if the flow field approaching, and out of, the grid was really 
steady then this auxiliary shock front will stabilize itself, de 
pending upon the effective grid solidity, at some suitable distance 
downstream of the grid 

The author believes that it is this auxiliary shock front which 
Glass and Patterson! have identified as a vortex in their position 
vs. time schlierens. The general agreement between their other 
results and ours” * and also a clear similarity of behavior of this 
auxiliary shock phenomena, viz. its appearance only at com 
paratively high incident shock strengths and its slow speed of 
advance, strengthens our belief that the said authors' were 
handicapped in analysis of this aspect of the experimental data 
by their experimental technique. A step-by-step record of the 
time history and development of the interaction by shadow 
graphs happens to be comparatively more revealing in this pres 
ent case . 

Behind the traveling reflected and transmitted shock waves, 
the flow field (in regions 3 to 5, Fig. 2) is riddled with diffracted 
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turbulent contact surface in the region 4. At higher incident MACH NUMBER M3 ae - in whi 
shock strengths the auxiliary shock fronts appear in region 4. ee 
This makes any complete theoretical analysis difficult to attempt. VS. as 
However, to extract as much information and/or inferences from { 
the experimental data as possible, some simplifying assumptions INCIDENT SHOCK STRE NGTH S alt 
have been resorted to. 7 
It is assumed that the reflected shock front very close to the Let 
upstream side of the grid becomes normal and plane—i.e., 50 7 T ——, madit 
parallel to the grid (not too far from reality, see Fig. 3a and also lo--- oe. .— 1... o- tions 
reference 3)—and that the flow in the regions 2 and 3 is steady, Pe ikes 
one-dimensional and isentropic. Extending the conventional coh } | spec! 
plane normal shock wave analysis to the reflected shock wave one 4 7 ad 
can then get - r in 3 
. ! 
M, = f(Usr, U,,...) (1) / dimen 
} / | re 
t 30 s a One al 
- ! The 
= P 
20K + 
° j | 
g | If one 
/ | shor k 
10 =. See ae a the c 
si relatic 
A 
For 
" 2 3 4. 5 € igh 
6 Wad ¢ I 
S—_—_ M iC h 
T oO 
Fic. 4 © Calculated points a. 
ment 
where 3 = u3/a3 = the Mach Number of the flow approaching one sl 
the grid in region 3 (called 17; in reference 1). “3 and a; represent proxit 
respectively the flow velocity and the sound speed in region 3 flow a 
U, = speed of the incident shock front Usr = speed of the Bef 
reflected shock front traveling upstream against the drift flow Seneill 
behind the incident shock wave. ae 
From the experimentally measured ambient conditions and of an 
Fic. 38a. Emergence of the auxiliary shock front through the shock speeds U, and U’sr, 3 can be calculated. It was observed? establ 
grid. © reflected shock wave; (_] auxiliary shock front; that, for the grid used, J; reaches an upper limit for S = 3.2 paper 
transmitted shock wave. S = 3.2. Delay time = 65 X 1076 and stays fairly constant for values of S higher than that (Fig tein 
os 4). Glass and Patterson! have also arrived at a similar conclu ences 
sion. 
It is important to point out that the auxiliary shock fronts Pay 
first appear only when the incident shock strength S is high est ; 
enough (depends on effective solidity of the grid). Since it is a 
observed that corresponding to these high values of S the Mach nae 
Number .43 assumes a constant value, obviously the grid is nt & 
choked.5 The grid holes act like crude nozzles. At their exit we 
the flow expands and develops a shock front similar to a terminal abe 
shock in a supersonic tunnel. The associated diamond-type pes 
wave pattern (see Fig. 3b) is similar in nature to the one ob — 
served, due to over-expansion, in the flow of a supersonic jet at 
the exit of a nozzle. 
Our experimental observation of this auxiliary shock front 
downstream of the grid is very well substantiated by the shadow : 
graphs of a supersonic ‘“‘Swiss Cheese’’ projectile fired on a 73-10 
ballistic range.® ay 
Let the solidity ratio of the grid be defined as ate 
livel 
s = 1— (A)/A De 
NACA 
where Ag is the projected open area of the grid exposed to the flow ‘Ke 
(in the present setup the grid is mounted perpendicular to the G 
flow direction), and A is the projected area of the grid boundary Lie 
exposed to the flow and is equivalent to the cross-sectional area dees 
of the shock tube. a Ut 
Fic. 3b. Emergence of the auxiliary shock front through the aan ey sation s te coneatinthy ties Seantional dngeee wipae = 
grid. [] auxiliary shock front. S = 3.2. Delay time = the grid obstructs the flow. Let us imagine that the solid and Fe 
1,375 X 10-6 sec. the open parts of the grid are so arranged that, at the section Th 
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in which the grid is mounted, they represent an idealized and 
symmetrical contraction and then enlargement of the cross- 
sectional area of the shock tube. The minimum cross section 
sav A,,) thus created is equivalent to Ao, the total open area of 
the grid 

Let us further assume that one-dimensional isentropic flow 
conditions are still valid. These are obviously crude assump- 
tions. However, the presence of an expansion through the grid 
makes the assumption of isentropic flow a bit less unpalatable, 


especially prior to the advent of the choking phenomena 


Under these assumptions, we have at our disposal the standard 
area ratio-- Mach Number relationship for an isentropic one 
dimensional channel flow.?. When the incident shock strength 
S(or M3) reaches a certain value, the flow reaches Mach Number 


One at the minimum cross section A,, of the shock tube. 
Therefore 


A.r/A = f(M (2) 


If one experimentally measures U, and U spr for some incident 
shock strength S for which the grid manifests choking, then 
the corresponding Mach Number 3 can be calculated from 
relation 1. Then by the use of relation 2, A.,/A and (since 
A Ay) the solidity ratio s of the grid can be determined. 


For the case in hand, the first indication of choking of the grid 
ippeared at incident shock strength S “ 3.0. The phenomena 
was clearly established at S = 3.2 (Fig. 3a). The corresponding 
Mach Number Mz; & 0.47 (Fig. 4). The calculated solidity 
ratio s = 0.29. The grid used had a solidity ratio s = 0.25 
(Fig. 1). Considering the bold assumptions involved, the agree- 
ment is surprisingly close. For a grid of known solidity ratio, 
one should, therefore, be able to predict within a tolerable ap- 
proximation, the value of Mach Number MM; of the approaching 


flow at which the grid will get choked. 


Before submitting this article for publication, the author re- 
ferred it to Glass and Patterson for their comments. They 
agreed that my recognition of the choked flow and the emergence 
of an auxiliary shock front downstream of the grid correctly 
establishes the true identity of the so-called ‘“‘vortex’’ in their 
paper.! The author first observed this choking phenomena 
during the summer of 1953 and the results are reported in refer- 


ences 2 and 3. 


Dr. Glass also disclosed that they are currently investigating 
the shock wave -—screen interaction by using schlieren and 
interferometer techniques. They furnished the author with 
some of their recent interferograms and schlieren photographs 
on the basis of which they too had recently recognized the true 
identity of the ‘“‘vortex’’. The author understands that Glass 
and Patterson expect to present the results of their recent inves- 
tigations of shock wave — screen interaction at the Fifth Inter 
national Meeting of the [AS—RAeS at Los Angeles. 
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FORUM 799 
Surge Behavior in a Three-Stage Compressor 


Carl E. Pearson 

Assistant Professor of Applied Mechanics, Harvard 
University, Cambridge, Mass 

June 8, 1955 


SUMMARY 


Experimental pressure and flow patterns corresponding to the onset of 
surge in a three-stage compressor have been studied under a variety of operat 
ing conditions I'wo very different types of surging, mild and deep, were 


encountered. The former always occurred in the presence of stall propaga 


tion and at an operating point which depended on the exit-chamber volume; 
its wave form was sinusoidal and its amplitude ranged from zero upwards 
The latter was of large amplitude, had a distorted wave torm, and corre 
sponded to complete flow blockage which was not overcome until the exit 
pressure had dropped to almost the inlet value. The onset of this deep surge 
was not affected by the exit volume, although its subsequent character was 


(1) CHARACTERISTICS OF COMPRESSOR 


oo COMPRESSOR consisted of a set of inlet guide vanes and a 
stator ring, followed by three stages, each of which contained 
a rotor and a stator. Inlet was atmospheric, and the discharge 
was diffused into a chamber of variable volume which in turn dis- 
charged through a calibrated exit throttle. The exit throttle was 
made up of a large number of small butterfly valves distributed 
over one end of the exit chamber 

A set of pressure and throttle characteristics is plotted in Fig. 
1 at a number of r.p.m. for the operating region near the peak of 
the pressure characteristics. It may be of incidental interest to 
remark that for flows lower than those corresponding to the pres- 
sure peaks, conventional flow pressure instrumentation was found 
to yield entirely erroneous operating points because of local flow 
irregularities. In the present compressor, the pressure drops 
abruptly as the flow is decreased beyond the pressure-peak point; 
the dotted portions of the pressure characteristics are unstable, 


and these portions cannot be retraced without some hysteresis 


(2) Tests with Exir CHAMBER VOLUME OF 46 CU.FI 


Two strain gage pressure pick-ups were placed 180° apart on 
opposite sides of the exit chamber; surge would then be evidenced 
by the coincidence in phase of the resulting tracings (Sanborn re- 
corder, trace speed 100 mm. per sec At all r.p.m., the traces 
were perfectly smooth until the flow was decreased beyond that 
corresponding to the peak of the pressure characteristic. When 
this point was passed by closing the throttle slightly, the charac- 
ter of the traces altered in a manner which depended on the r.p.m 
Figs. 2 and 3 show typical traces for r.p.m. values less than and 
greater than 4,000 r.p.m., respectively 

The surging of Fig. 2 is called mi/d and that of Fig. 3 is called 
deep. The former could be of any amplitude down to zero, 
whereas the total pressure fluctuation of the latter was of the 
same order as the difference between peak pressure and atmos- 
pheric pressure. (For example, Fig. 3 was taken at 5,500 r.p.m., 
and the double amplitude of the trace was 0.11 atmospheres 
The pressure scale in Fig. 2 is magnified by a factor of four as 
compared to Fig. 3.) 

At the very lowest r.p.m., mild surging did not begin imme- 
diately the pressure peak points was passed; as the throttle was 
closed further, mild surge began in an intermittent fashion, and 
became steady and of increasing amplitude with still further 
closure. Finally, with the throttle almost completely closed, 
deep surge ensued. As the r.p.m. increased (but remained be 
low 4,000), mild surge began sooner with throttle closure; at 
4,000 r.p.m., it began immediately the pressure peak point was 
passed 

The presence of stall propagation! at various flows was next 
examined by placing two (uncompensated) hot wires 90° apart 
behind the first rotor, and again taking simultaneous tracings 
No stall propagation was observed at any r.p.m. until the flow was 
reduced beyond the peak pressure point. For r.p.m. below 4,000, 
definite stall propagation then occurred for all lower flows (and it 
is of course this alteration in the flow character which causes the 
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Fic. 1. Portion of pressure and throttle characteristics of three-stage compressor. 
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Fic. 2. Mild surging Fic. 3. Deep surging 
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;dden drop in total pressure A typical trace is shown in Fig 


rhe time separation between the two traces corresponded to 


sual stall propagation speeds rhe frequency decreased as the 
isu I ] 


was further closed For r.p.m. above, 4,000, the deep 


throttle 


masked the stall propagation known from previous work 


surge 


to be present during a portion of each surge cycle 


In deep surge, that portion of the curve corresponding to a 


ire decrease was found to be parallel to the throttle dis 


pressl 


harge line of the exit chamber—i.e., during this portion of the 


surge CVC le, the compressor was « ompletely blocked rhis result 
was checked by taking simultaneous tracings of the exit chamber 
voltage across a hot wire placed ahead of the 
Fig. 5 The 


part of the hot-wire trace corresponds to about 10 per cent of the 


pressure and the d.c 


first rotor. A typical tracing is shown in lowest 


velocity represented by the upper part of the trace (the vertical 
scale is not linear, the wire being most responsive at low veloci 
ties 

By setting the operating point on the verge ol deep surge, and 
then momentarily slightly closing and reopening the throttle, as 
In no case 
had 


occurred, it could not be overcome until the exit pressure attaine 1 


little as a single surge pattern resulted (see Fig. 6 


could less than this pattern be obtained—i.e., once blockage 
i level near atmospheric It was found that the pressure fall-off 
time in deep surge was always about the same, corresponding to 
the calculated time required to discharge the exit chamber down 
toa level near atmospheric (the precise level affects the discharge 
time only logarithmically It will be noticed for example in Fig 
7, which shows the effect of running into deep surge along a con 
stant throttle line, that although the period of the surging is much 
The 


7 is probably to be found in 


longer than in Fig. 8, the fall-off time is about the same 
reason for the longer period of Fig 
the more delicate surge point setting available by use of the cou 
stant throttle line; it takes a long time for the exit pressure to 


build up to the blocking value 


3) TESTS WITH Exit CHAMBER VOLUME OF 12 CU.F1 


Phe behavior at speeds below 4,000 r.p.m. was similar to that 


described in Section 2 as regards stall propagation and surging, 
except that mild surging was weaker and required more throttle 
For above 4,000, the initial deep 


closure for inception r.p.m 


surge was as before, except that it now did not continue. A 


typical trace is shown in Fig. 8 (same operating point as Fig. 3 
and 
The 


subsequent average pressure level was lower than with the larger 


Again, the fall-off times were roughly the same for all r.p.m., 
corresponded to the discharge time for the new volume 
volume. After the initial deep surge cycle, mild surging of various 


unplitudes was found superposed on stall propagation 


DISCUSSION 


If the air were perfectly incompressible, then the same propor 
tional increase in both flow and r.p.m. would not alter the interior 
flow angles, and the ratio of r.p.m. to flow corresponding to the 
pressure fall-off point would be the same at all r.p.m For the 
ibove compressor, this result is indeed approximately correct 
Compressibility does however have the effect that (because of 
general design considerations) the last stages are closer to stall 
at the higher r.p.m. than the lower, when the first stage stalls 


rhis means that stall should be more severe at the higher r.p.m 


Phere is an associated effect involving the pumping action of 
the last stages. It is known! that the same rotor will permit a 
variety of stall propagation patterns to exist at a particular flow 
nd r.p.m. point, the total area stalled increasing as the pumping 
iction of the latter stage decreases. In the present machine, the 


less-compressed flow at the lower r.p.m. has a relatively greater 


ixial velocity in the latter stages, and so the pumping action of 
these stages is affected by first stage stall less than at the higher 
rp.m., thus favoring less stalling in the first stage 

The combination of these two related effects would seem to 


count for the difference between the character of the pressure 
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fall-off encountered in the twor.p.m. ranges. It may incidentally 
be pointed out that the slope of the pressure characteristic at that 
point at which deep surge commences is so steep that this term 
dominates all others in the criteria of reference 2 and in fact en- 
sures fulfillment of the criteria. The question now arises as to 
why the deep surging continues with large volumes but not with 
small. Perhaps the explanation lies in the fact that with the 
larger volume, the pressure fall-off time and the beginning of the 
rise time are sufliciently long for the compressor to thoroughly 
purge itself of stalled flow, whereas this is not the case for the 
shorter times associated with the smaller volume so that the com 
pressor in the latter case falls naturally into a stall propagation 
pattern. The effect would seem to be similar to that of hysteresis 
in a compressor, in which the flow may have a different character 
depending on whether it is approached from the high or low flow 


sides. 


Mild surge in the compressor always takes place in the presence 
of stall propagation, but stall propagation may exist without surge 
Presumably, the locally steep slope of the pressure characteristic 
required by reference 2 for surge is obtained by the possibility of 
alteration in the size of the stalled regions propagating around the 
wheel. Its incidence and character is very much affected by the 
exit chamber volume, as predicted by reference 2 

Insofar as the results of these tests may be generalized to apply 
to jet engine compressors, the conclusions of reference 2 may be 
extended by the following remarks: 


(a) Other than the position of the operating point, there seems 
to be no overall effect-——such as the presence of pressure fluctua 
tions which warns of the imminence of the pressure fall-off 
point. (This statement may have to be modified for those com 
pressors which stall-propagate before the fall-off point is reached, 
where the character of the stall propagation might provide such 


a clue.) 


(b) During acceleration, the tendency to surge is greatest in the 
vicinity of the knee of the surge line; this corresponds to a 
moderate value of r.p.m. which conceivably could fall either side 
of the r.p.m. dividing line described above, depending on the en 
gine design. It is then possible that the designer may exercise 
some choice as to the type of pressure fall-off —i.e., as to whether 
the pressure is to fall to a steady nonoscillatory lower value well 
above atmospheric or almost to atmospheric with at least initial 


oscillations. 


(c) It also appears that there is a critical value of combustion 
chamber volume, such that (assuming the r.p.m. sufficiently high ) 
for smaller volumes only one surge cycle will occur whereas for 
larger volumes continued surging takes place. The latter solu 
tion may be mechanically and thermally destructive, but does 
have the advantage of maintaining a higher average pressure 


(d) If surge occurs at fairly high r.p.m., then the compressor 
remains totally blocked until its exit pressure has fallen to a cer 
tain level, and the only way in which recovery can be speeded is 
to increase the discharge rate of the exit air. The important 
fact seems to be that the blockage process can not be interrupted, 
as, for example, by using an emergency discharge nozzle to open 
automatically when the exit pressure begins to fall. 


Because of the way in which the slope of the characteristic 
dominates the surge criteria of reference 2, it follows that the only 
way of preventing surge of this stall-blockage type is to inhibit the 
stall itself by some of the methods there discussed; alteration in 
burner volume or rate of fuel-burning, for example, would be ex- 
pected to have little effect 
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Supplementary Remark on ‘‘Airplane 
Stopping Distance’’’ 


B. Saelman 
Design Weight Engineer, Lockheed Aircraft Corp., Burbank, Calif 
June 9, 1955 


I REFERENCE | a formula was given for airplane stopping dis 
tance. It should be noted that the basic differential equation 
for the landing ground run is integrable when the following devices, 
not mentioned in the above reference, are taken into account 

(1) Jet reverse thrust, usually considered constant, but may 
be considered proportional to a power of ground speed 

(2) Downhill or uphill component of gravity, constant 

(3) A variable coefficient of friction expressible as a poly- 
nomial in ground speed. 

(4) Drag chute 

For example, if 


“w= po — C2 (1 


and the jet reverse thrust is constant, the differential equation 


of motion reduces to 


(A — By)y + Dy? + Ey? — Fy + G = 0 (2 
where 
y = Ff 
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l l 
and 
Co, = chute drag coefficient 
Ss, = drag chute area 
6 = runway slope 
J = jet reverse thrust 
KiJ =0 
K,J = 0 


where the nomenclature is the same as in reference 1. 


Eq. (2) can be expressed 
By—A 


ly = dl (3) 
Dy? + Ey?— Fy+G ‘ 


The coefficient of dy is a rational fraction and can be integrated 
by the method of reference 2, assuming that it can be resolved 
into a series of partial fractions 

From the brake characteristics, a curve of coefficient of friction 
versus velocity can be obtained. By substituting sets of points 
from the characteristic brake curve into the general polynomial 
relationship, an analytic expression for « versus .¢ can be obtained 
This can then be used in the above analysis for stopping distance 

The method provides a refinement only if the actual variation 
in the coefficient of friction is described by Eq. (1) or a poly 
nomial in x. When the coefficient of friction is described by 
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different functions over the various segments of the run, an 
integration by segments may be performed 

Possible additional refinements would include the effects of tail 
loads, noncoincidence of lift and gravity forces, induced drag re 
sulting from lift and ground effects, and rolling drag resistance of 


nose W heel 
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Measurements of the Effect of Two- 
Dimensional and Three-Dimensional 
Roughness Elements on Boundary-Layer 
Transition* 

P. S. Klebano#, G. B. Schubauer, and K. D. Tidstrom 

Fluid Mechanics Section, National Bureau of Standards, 
Washington, D.C 

June 15,1955 


- HIs stuby of the effect of roughness on transition, H. L 
Dryden! found, on the basis of available data, that the effect 
of a two-dimensional roughness clement such as a “trip wire”’ 
could be represented reasonably well in terms of a functional re 
lation between Re, and k/6,*, where Re, is the Reynolds Number 
of transition based on distance from the leading edge, & is the 
height of the roughness element, and é,* is the boundary-layer 
displacement thickness at the position of the element. At his 
suggestion some additional data were obtained, primarily to 
extend the range to higher values of k/6,*, during the course of an 
investigation of transition on a flat plate conducted at the 
National Bureau of Standards.? After the results on the two 
dimensional roughness elements were obtained, it appeared to be 
desirable to see whether a row of three-dimensional roughness 
elements would behave in the same way 

The roughness elements were ‘attached to the smooth flat 
plate with sharpened leading edge on which transition studies 
were conducted. The plate spanned the test section of the 
Bureau's 4!/.-ft. wind tunnel in which the turbulence level was 
sufficiently low for the Reynolds Number for the start of transi 
tion without roughness, Reo, to be 2.5 K 10° at a wind speed of 
80 ft. per sec. A small favorable pressure gradient existed over 
the first 1'/, ft., but from there on the pressure was constant to 
within +0.5 per cent of the dynamic pressure 

The two-dimensional roughness elements were cylindrical rods 
which were cemented individually to the surface 2 ft. from the 
leading edge. They spanned the plate for a distance of 1'/, ft 
on each side of the centerline and ranged in size from 0.033 in 
to 0.1875 in. The three-dimensional elements were spheres 


*Sponsored by the NACA 
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cemented to the surface in a row perpendicular to the flow at 
2 ft. and also 1 ft. from the leading edge with various spacings 
Spheres !/32-, '/is-, and '/,-in. in diameter were studied Fig. | 
shows the experimental arrangement, which in this case consisted 
of the '/i.-in. spheres with a spacing twice the diameter and ex 
tending for a distance of 6 in. on each side of the centerline 
Extension of this distance had no effect on the results. The 
hot-wire anemometer used to detect transition is shown in posi 
tion for traversing along the centerline downstream from the 
spheres 

It is well known that transition need not occur at the rough 
ness element. If the speed is low enough or the height of the 
roughness is small enough, transition will take place some 
distance downstream, beginning with the appearance of small 
patches of turbulence which develop in size and number over a 
transition region to form the fully turbulent layer, as described 
in reference 2. As the speed or the height is increased, the re 
gion moves upstream and becomes shorter until it finally reaches 
a limiting position at the roughness. It is then said to be at 
tached. As shown in reference 1, the limiting position defines 
limiting lines in Fig. 2 expressed by 


Re, = 2.96 (xx/k)? (k/6,*)? 


where x; is the distance of the roughness element from the lead 
ing edge. These limiting lines apply to both two- and three 
dimensional elements, but are shown in Fig. 2 by the short 
dashed curves only for the three-dimensional results. All of the 
plotted results pertain to the unattached condition, and values 
of Re, refer to the beginning of transition. Depending on the 
size and position of the roughness clements, the permissible speeds 
ranged from 6 to 80 ft. per sec 

Dryden found that Re,/ Rep correlated data obtained in different 
wind tunnels better than did Re, alone, account thus being taken 
of variations in turbulence from tunnel to tunnel. However, 
Re, is used here, since the present results pertain to only one 
wind tunnel and since in reference 3 it is shown that the effect of 
stream turbulence becomes negligible for values of k/6,.* > 0.6 
Fig. 2 shows in marked contrast the difference between the 
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effects of two- and three-dimenusional roughness. Whereas the 
two-dimensional results indicate a relationship between Re, and 
k/6;.* in agreement with Dryden’s findings, the three-dimensional 
results fail to satisfy such a relationship. The three-dimensional 
than the 


* 


roughness has less effect on transition for a given k/6, 
two-dimensional roughness, but the effect is of a more critical 
nature. This is seen from the very narrow range of k/6;* cov 
ered by a given size, indicating a rapid change in transition 
Reynolds Number for small changes in wind speed. Fig. 3 
shows an attempt to represent the effect of the spheres in terms 
of a roughness Reynolds Number, l.k/v, where LU; is the ve 
locity at the height of the sphere. The scatter is appreciable 
and the nearly vertical array of points shows the sensitivity of 
Re, to Uyk/v once the effect of the roughness begins to be felt 
This means that the roughness is essentially without effect until 
a value of L’.k/v is reached which is likely to bring transition to 
the point of attachment. The following critical values for at 
tachment appear therefore to be the practically important quan 
tities: 
(U,k/v), 

16-in. spheres at 2 ft., spacing Sk 573 

16-in. spheres at 2 ft., spacing 4k 190) 

16-in. spheres at 2 ft., spacing 2k 518 

'8-in. spheres at 2 ft., spacing 2k 760 

8-in. spheres at 1 ft., spacing 2k 890 

32-in. spheres at 2 ft., spacing 4k 542 


32-in. spheres at 1 ft , spacing 4k 578 


Excluding the value of 890, the critical values are in reasonable 
agreement and the average is 577. The corresponding value 
found for the two-dimensional roughness lies between 200 and 
300. It may be of interest to note that the ratio of approxi 
mately two corresponds to the relative perturbations of spheres 
and cylinders in potential flow. 

The value of 890 is found for the curve at the extreme right 
of Fig. 2. Here the curve follows a pseudo limiting line, corre 
sponding to a fixed distance downstream from the roughness, 
and then finally drops to the real limit. This behavior is appat 
ertly associated with the stability of the layer, and it is note 
worthy that the intersection with the real limiting line does not 
occur until this line itself rises above what appears to be an 


asymptotic limit of Re, ~ 90,000. The corresponding Rey- 
nolds Number based on displacement thickness is 520, which js 
about the critical value found for a Blasius distribution for sta- 
bility against small perturbations. Apparently self-sustaining 
turbulence will not develop below this limit no matter how large 
the initial disturbance. One may question why this limit should 
apply when obviously the mean velocity is altered just back of 
the roughness. It is observed that when transition occurs some 
distance from a roughness element, the alterations in the velocity 
profile have disappeared, leaving, however, perturbations that 
would not have been present in the absence of the roughness 
It appears, therefore, that the primary effect of the roughness 
is to increase the magnitude of the perturbations, probably 
through the unstable character of the flow behind the roughness 
This argument does not hold when transition occurs at the ele 
ment, but in the case of a three-dimensional element, the develop- 
ment of a wedge of turbulence starting from the element evi- 
dently is not possible when the surrounding flow is below the 
stability limit.2 This is a point which may possibly deserve 
further investigation in view of the possibility of obtaining much 
higher stability limits by means of distributed suction 
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Errata—‘‘On Buckling of a Clamped Rhombic 
Thin Plate in Shear’’! 


M. Hasegawa 
Shimane University, Matsue, Shimaneken, Japan 
June 22, 1955 


I REFERENCE I, the following errata have been found 
The sentence immediately preceding Eq. (1) should read: 


“If w be the transverse displacement of a point on the middle 
plane, the total energy Ll’ of the clamped rhombic plate is ex- 
pressed by the following: (see Fig. 1.)"’ 

Eq. (1) should read: 


D ve i O7u' 
2U = — (A,w)? $ cosa Ayu 
sin? a. ‘ . | Oxoy 


ow \?] 2D f* "at Ow Ow 
t cos? a dxdy + dxdy 
Oxoy f ee 3 ; , OX OY 


Fig. 1 should be inverted. 
Table 1 should read: 


60° 

98.99 
The last paragraph should read: ‘‘The true value of Xo cal 
culated by the author for square plate is 36.13, and B. Budiansky 
and R. W. Connor found that 36.12 S A» S 36.49. By compar- 
ing these values, the values in the table seem to be of good ac 


curacy.” 
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